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APPENDIX  B. 


1  .  INTRODUCTION 


For  several  years,  A.R.A.P.,  Inc.  has  been  developing  a  computer  model 
for  determining  the  detailed  low-level  atmospheric  distributions  of  velocity, 
temperature,  moisture,  refractive  index,  and  the  turbulent  variations  of  these 
quantities  for  marine  environments.  The  three  physical  processes  most 
critical  for  determining  the  atmospheric  marine  boundary  layer  are  turbulent 
transport,  thermal  radiation,  and  change  of  phase  of  atmospheric  water. 
Reference  1  provides  a  review  of  our  understanding  of  these  processes  and  a 
review  of  some  of  the  sample  calculations  which  successfully  illustrate 
features  expected  in  the  atmospheric  marine  boundary  layer.  The  most  recent 
status  report  on  this  model  is  given  in  Reference  2.  Details  of  the 
foundation  of  the  model,  yearly  developments  and  a  number  of  exemplary 
simulations  are  given  in  References  3-19. 

Our  efforts  over  the  past  18  months  have  again  been  divided  between 
simulations  using  the  existing  model  and  model  developments.  A  major  part  of 
our  model  simulations  have  been  associated  with  the  particular  scenarios 
chosen  by  Calspan  for  use  in  their  fog  model  study.  These  simulations,  along 
with  the  model  modifications  made  as  a  result  of  these  calculations,  are 
discussed  in  Chapter  2.  Another  important  segment  of  model  simulations  has 
been  devoted  to  further  understanding  the  process  of  wave  breaking  in  stably 
stratified  regions  of  the  atmosphere.  The  goal  of  these  simulations  discussed 
in  Chapter  3  is  to  provide  an  adequate  representation  of  this  wave-turbulent 
interaction  for  one-dimensional  models  of  the  atmospheric  boundary  layer. 

The  two  principal  model  developments  during  this  time  period  are  improved 
modeling  of  the  cloud  microphysics,  and  the  incorporation  of  the  possibility 
of  an  anisotropic  scale  into  our  second-order  closure  model.  These  two 
developments  are  presented  in  Appendix  A  and  B,  which  are  written  for 
submission  as  separate  journal  articles.  Preliminary  explorations  of  two 
major  model  developments  are  presented  in  Chapters  4  and  5.  These  are  the 
development  of  a  three-dimensional  model,  and  cumulus  parameterization. 


2.  PARTICIPATION  IN  THE  CALSPAN  FOG  MODEL  EVALUATION  STUDY 


Introduction 

As  part  of  an  assessment  of  the  Naval  Air  System  Command's  Marine  Fog 
Investigation  program.  Cal  span  organized  an  evaluation  study  of  fog  models 
based  on  comp?irison  of  model  simulations  for  six  observational  cases  chosen  by 
Calspan.  A.R.A.P.  participated  in  this  evaluation  study  by  first  responding 
with  the  "blind"  model  simulations  requested  by  Calspan,  and  second  by  a  rerun 
of  each  of  these  cases  after  the  preliminary  results  were  made  public  at  the 
workshop  at  NEPRF  in  May  1982. 

The  initial  model  results  were  rather  disappointing  for  us,  with  only 
three  of  the  cases  showing  reasonable  correspondence  with  data.  Our  prime 
purpose  here  is  to  report  on  subsequent  model  developments  and  the  information 
gained  from  model  reruns  of  these  fog  scenarios. 

Brief  Review  of  the  A.R.A.P.  Model  Operation  for  this  Fog  Study 

We  used  the  one-dimensional  version  of  our  Reynolds-stress  transport 

closure  model,  in  which  differential  equations  are  solved  for  all  first  and 

second-order  moments  of  the  dynamic  and  thermodynamic  variables.  The 

1  12 

equations  are  described  in  detail  by  Lewellen  ,  and  by  Oliver  et  al.,  .  The 
differential  equations  are  solved  numerically  using  finite-difference 
techniques  and  a  self-adjusting  grid  system,  i.e.,  the  numerical  mesh  changes 
as  the  solution  changes  to  maintain  good  resolution  of  sharp  features  such  as 
the  capping  inversion. 

The  six  Calspan  fog  cases  were  simulated  by  marching  forward  with  time  as 
the  independent  variable,  and  converting  the  spatial  surface  temperature 
variation  into  a  temporal  variation  using  the  mean  velocity  in  the  boundary 
layer.  We  recognize  that  this  does  not  yield  a  completely  consistent 
simulation  of  the  data  which  represents  vertical  distributions  of  temperature 
and  humidity  at  different  horizontal  positions  at  different  times.  However, 
it  probably  represents  as  consistent  a  simulation  as  is  possible  using  data 
which  are  insufficient  to  specify  the  effects  of  horizontal  advection.  Our 
simulation  thus  necessarily  neglects  contributions  from  horizontal  advection, 
which  are  very  likely  to  be  a  significant  source  of  error. 
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The  atmospheric  radiation  scheme  treats  terrestrial  long  wave  and  solar 
shortwave  radiation  separately.  The  long  wave  component  is  calculated  using 
the  usual  convolution  integrals  over  the  emission  and  transmissivity  profiles. 
The  transmissivities  are  parameterized  as  sums  of  exponential  functions  of  the 
liquid  water,  water  vapor,  and  carbon  dioxide  path  integrals.  The  emission 
functions  are  taken  to  be  the  black  body  emission  at  the  local  absolute 
temperature  of  the  air. 

Shortwave  radiation  is  calculated  using  the  "two-stream"  model  for  upward 
and  downward  flux  components.  The  scheme  includes  a  direct  absorption 
coefficient,  and  also  a  scattering  effect  when  liquid  water  is  present. 
Details  of  the  radiation  schemes  can  be  found  in  Lewellen  et  al.,^**. 

Modifications  Arising  After  the  NEPRF  Meeting 

As  a  result  of  the  comparison  between  initial  model  results  and 
observations  shown  in  Figures  1  to  12  and  presented  at  the  meeting  at  NEPRF  in 
May  1982,  we  made  an  extensive  review  of  our  model.  The  program  was  modified 
to  correct  several  possible  problems,  and  all  six  cases  rerun  with  the 
improved  model.  In  addition  to  correcting  some  input  conditions,  the 
following  program  changes  were  made: 

(i)  The  full  temperature  profiles  through  the  troposphere  and  higher  have 
been  incorporated  into  the  downward  long  wave  radiative  flux 

calculation.  This  eliminates  the  need  to  estimate  downward  fluxes  at 
the  top  of  the  boundary  layer  from  the  profile;  the  flux  is  now 
explicitly  calculated  using  the  same  radiative  transfer  model  as  the 
boundary  layer  program. 

(ii)  The  reflection  of  long  wave  radiative  flux  at  the  ground  has  been  set 
to  zero.  In  our  original  runs,  there  was  a  reflection  at  the  surface 
of  10%  of  the  downward  radiation.  This  produced  unreasonable  heating 
close  to  the  ground  in  the  stable  cases  2  and  5,  and  seems  to  be  due  to 
the  inconsistency  between  assuming  black  body  emission  properties  and 
non-zero  reflectivity. 
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(m)  A  numerical  finite-differencing  problem  at  the  top  of  the  cloud  layer 
was  pointed  out  by  Steve  Burk  where  the  radiative  fluxes  are  virtually 
discontinuous.  Our  use  of  a  central  difference  resulted  in  part  of  the 
large  cooling  being  applied  to  the  dry  layer  immediately  above  the 
cloud.  This  produces  a  spurious  entrainment  of  the  dry  air  into  the 
boundary  layer.  This  problem  was  corrected  by  using  one-sided 
differences,  and  care  was  taken  to  ensure  that  the  total  flux  budget 
was  preserved. 

(iv)  The  approximate  formula  for  the  radiative  convolution  integral  as 
described  in  Lewellen  et.  al,^^  was  replaced  by  the  full  integral 

expression;  this  is  more  expensive  to  compute  but  should  be  more 
accurate. 
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Model  Results 


Case  1 


This  case  corresponds  to  observations  during  a  cruise  Northwest  of  San 
Nicolas  Island  on  22  May  1978.  The  simulation  covers  4  hours  from  0500  POT  to 
0900  PDT.  Air  temperature  was  colder  than  the  sea  surface  temperature.  This 
case  was  run  with  two  different  initial  conditions,  first  with  a  thin  cloud 
layer,  then  with  no  initial  cloud.  In  both  cases  the  cloud  base  lowers  and 
the  cloud  becomes  thicker.  After  four  hours,  the  results  are  not  very 
dependent  on  the  initial  conditions.  No  fog  was  observed  or  predicted.  The 
new  model  results  for  temperature.  Figures  13  and  14,  are  slightly  closer  than 
the  old  results,  but  both  runs  were  in  reasonable  agreement  with  the 
observations.  The  major  difference  in  our  most  recent  run  and  the  initial  run 
was  a  correction  to  the  input  solar  radiation  which  was  erroneously  set  at 
only  4%  of  its  correct  value. 

Case  2 

This  case  calls  for  the  simulation  of  the  shallow  advection  fog  formed 
over  the  cold  sea  surface  during  a  3  hour  evening  period  (1700  to  2000  EOT) 
off  the  coast  of  Nova  Scotia  on  2  August  1975.  The  new  run  for  this  case. 
Figure  15,  is  much  closer  to  the  observations  than  our  original  result.  The 
ground  reflection  of  long  wave  radiation  was  the  cause  of  the  low-level 
heating  in  our  earlier  integration  and  gave  temperatures  roughly  3°C  too  warm. 
Without  the  reflection,  our  predicted  temperatures  are  very  close  to  the 
observations  below  30  m  where  the  measurements  were  made.  We  still  do  not 
predict  fog,  but  our  maximum  relative  humidity  is  now  about  96%  in  the  lowest 
few  meters.  We  believe  the  relatively  large  horizontal  gradients  in  sea 
surface  temperature  observed  at  this  time  (Figure  16)  play  a  key  role  in  this 
type  fog.  As  mentioned  previously,  the  simulation  of  this  as  an  unsteady, 
horizontally  homogeneous  flow  should  not  be  expected  to  accurately  simulate 
the  role  of  horizontal  advection  of  temperature  and  humidity.  This  could  only 
be  the  case  if  vertical  gradients  in  wind  velocity  were  unimportant. 
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Figure  2.13  Case  1,  Cloud,  T  =  4  hours 
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Figure  2.14  Case  1,  Cloud,  T  =  4  hours 
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Figure  2.15  Case  2,  T  =  3  hours 
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Figure  2.16  Isopleths  of  sea  surface  temperature  (°C)  and  regions  of  fog 
demonstrating  the  spatial  inhomogeneities  present  during  the 
case  2  simulation. 


Case  3 


This  case  is  a  simulation  of  20  hours  observed  between  San  Nicolas  Island 
and  San  Diego  on  7  October  1976.  A  100  m  thick  stratus  cloud,  beneath  a  250  m 
inversion,  dissipates  during  the  day  and  redevelops  during  the  early  evening 
to  produce  fog  under  conditions  of  strong  subsidence  above  the  inversion.  Our 
original  run  for  this  case  was  grossly  in  error  because  the  shortwave  solar 
radiation  was  set,  inadvertently,  at  4%  of  its  correct  value.  Since  this  case 
runs  throughout  the  day,  solar  heating  is  quite  important,  and  the  differences 
in  the  new  run  are  mainly  due  to  this  change.  For  the  three  verification 
times,  we  see  that  the  model  predicts  a  shrinking  boundary  layer  in  response 
to  the  applied  subsidence,  and  also  significant  heating  during  the  daytime. 

The  run  is  initialized  with  a  cloud  about  100  m  thick  and  base  at  175  m. 
This  cloud  actually  lowers  to  form  a  fog  around  dawn,  but  is  then  evaporated 
by  the  solar  heating  over  the  next  2-3  hours,  so  that  at  the  first 
verification  time  of  10  hours.  Figure  17,  we  have  no  cloud.  At  this  time  our 
prediction  is  somewhat  cold  in  the  mixed  layer,  and  the  mixed  layer  depth  is 
too  shallow.  However,  at  the  next  verification  time  of  14  hours.  Figure  18, 
the  prediction  is  in  almost  exact  agreement  with  the  observations.  It  seems 

that  the  subsidence  rate  was  not  uniform  over  the  first  14  hours  as  assumed  in 
the  model  run. 

The  prediction  is  also  in  error  at  the  final  time  of  20  hours.  Figure  19, 
where  the  observations  show  a  fog  layer  up  to  250  m.  Our  prediction  gives  a 
70  m  boundary  layer  with  a  97.5%  relative  humidity.  There  are  two  possible 
causes  of  the  discrepancy.  Firstly,  the  assumption  of  constant  subsidence 
rate  may  be  incorrect;  it  may  well  be  reduced  in  the  later  part  of  the  run 
since  there  is  no  data  above  the  inversion  to  indicate  a  continued  downward 
trend.  Indeed,  to  entrain  over  150  m  over  6  hours,  i.e.,  from  100  m  to  250  m 
between  14  hours  and  20  hours  in  the  face  of  a  0.3  cm  s  ^  subsidence  at  that 
height,  seems  a  formidable  entrainment  rate.  The  second  possibility  is  that 
the  predicted  entrainment  rates  are  too  small;  this  is  quite  possible  also 
since  we  are  aware  of  deficiencies  in  the  model  entrainment  in  the  free 
convection  regimes.  The  discrepancy  is  most  probably  a  combination  of  the  two 
effects,  but  without  an  observational  measurement  of  one  of  them  it  is 
impossible  to  determine  their  relative  importance. 
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Figure  2.17  Case  3,  T  =  10  hours 
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Figure  2.18  Case  3.  T  =  14  hours 
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Figure  2.19  case  3.  T  »  20  hours 
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Case  4 


This  overnight  case,  2000  POT,  14  July  to  0530  POT,  15  July,  1973, 
simulates  a  cruise  80  km  offshore  Northwest  of  Pt.  Conception.  The  coastal 
stratus  thickens  from  2100  m  to  ''  350  m  depth  under  an  inversion  at  500  m. 
Although  our  original  integration  for  this  case  was  quite  close  to  the 
observations,  it  must  be  noted  that  we  needed  to  apply  a  significant 
subsidence  to  prevent  the  mixed  layer  from  growing.  This  was  a  result  of 
spurious  entrainment  arising  from  the  finite  difference  approximation  for  the 
radiative  flux  divergence  as  mentioned  earlier.  The  new  run  has  no 
subsidence,  and  does  not  give  boundary  layer  growth.  The  new  run.  Figure  20, 
is  even  closer  to  the  observations;  in  fact  the  agreement  is  almost  exact. 

Case  5 

This  is  another  shallow  advection  fog  formed  over  cold  water  ~  80  km 
offshore  Southeast  of  Nova  Scotia  on  5  August  1975.  A  3  hour  time  period  is 
simulated  starting  1  hour  after  sunrise.  The  fog  shows  a  dramatic  increase  in 
depth  when  it  flows  over  substantially  warmer  water.  As  in  Case  2,  our 
original  run  for  this  case  was  dominated  by  the  reflection  of  long  wave 
radiation  at  the  ground.  With  the  reflection  removed,  our  predictions  are  now 
much  closer  to  the  observations.  The  predictions  are  still  somewhat  too  warm. 
Figures  20-22,  with  the  result  that  our  relative  humidities  are  too  low,  and 
no  fog  is  formed. 

We  note  that  the  fog  in  this  case  is  formed  during  the  first  hour,  before 
the  surface  temperature  has  changed.  Inspection  of  the  profiles  at  t=l  hour 
suggests  that  the  water  vapor  has  been  brought  down  from  aloft  where  it  was 
initially  higher  than  the  surface  by  some  20%.  The  model  fails  to  predict 
this  transport.  We  initialize  the  run  with  super-equilibrium  turbulence  in  a 
stable  layer,  which  implies  very  little  mixing  away  from  the  surface; 
however,  the  problem  is  not  solved  by  specifying  higher  initial  turbulence 
levels.  The  profiles  at  t=l  hour  show  the  curious  phenomenon  of  humidity 
being  mixed  down  from  aloft  but  not  temperature.  When  the  model  turbulence 
levels  are  increased,  both  quantities  are  transported  downward  with  very 
little  change  in  relative  humidity,  and  no  fog  formation.  In  view  of  the  lack 
of  fog  in  this  run',  it  is  not  surprising  that  our  predicted  temperatures  are 
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Figure  2.20  Case  4,  T  =  10  hours 
- INITIAL 


8 .  VERIF 

0 - model 


27 


Figure  2.21  Case  5,  T  =  1  hour 
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Figure  2.22  Case  5,  T  =  2  hours 
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too  high,  since  we  do  not  nave  the  fog-top  long  wave  radiation  cooling  which 
must  have  been  present  in  the  observations. 

Case  6 

This  case  simulates  3  hours  just  before  sunset  during  which  stratus  was 
observed  to  thicken  and  produce  fog  in  a  marine  layer  capped  at  a  height  of 
200  m  in  the  Farallon  Islands  on  29  August  1972.  Although  our  earlier  runs 
for  this  case  were  in  good  agreement  with  the  observations,  the  new  runs  are 
significantly  different.  The  main  cause  of  the  difference  is  the  corrected 
value  for  the  solar  constant.  The  sunset  is  at  t=2.5  hours,  while  the 
verification  time  is  t=3  hours.  It  appears  that  our  formulation  for  the 
shortwave  scattering  produces  too  much  heating  in  this  case,  because  the 
heating  exceeds  the  long  wave  cooling  in  the  early  stages  and  quickly 
evaporates  the  cloud.  A  possible  cause  is  the  failure  to  account  for  the  low 
solar  angle,  which  reduces  the  shortwave  absorption  as  the  sun  approaches  the 
horizon,  according  to  Stephens  (Ref.  42,  Figure  3).  The  effect  depends  on  the 
state  of  the  atmosphere  above  the  cloud,  and  has  not  been  examined  in  any 
detail,  but  the  remedy  would  clearly  involve  making  the  absorption 
coefficients  dependent  on  the  solar  angle. 

The  temperature  profiles  at  t=3  hours,  Figures  23  and  24,  show  that  the 
prediction  is  too  warm;  this  is  as  expected  since  the  cloud  has  evaporated 
and  the  cooling  mechanism  has  been  lost.  A  run  has  been  made  without  the 
shortwave  radiation,  and  this  shows  the  cloud  descending  but  does  not  reach 
the  ground  in  3  hours.  Examination  of  the  details  of  the  integration  reveals 
some  sensitivity  to  the  initial  development  in  this  case.  At  t=0,  we  have  a 
region  of  high  humidity  in  the  cloud,  and  lower  humidity  below.  As  the  cloud 
top  cools  (without  the  sun),  the  entire  boundary  layer  mixes,  and  reduces  the 
humidity  in  the  cloud.  The  timing  of  this  mixing  event  is  important,  in  that 
if  it  is  delayed  then  there  is  more  cooling  of  the  layer  prior  to  the  cloud 
humidity  reduction.  In  our  model  runs,  the  turbulence  builds  up  and  mixes  the 
layer  in  30  minutes  or  less,  even  if  we  begin  with  very  low  turbulence  levels. 
At  this  stage,  the  cloud  thickness  is  drastically  reduced,  and  the  cloud 
becomes  optically  thin  since  it  has  not  cooled  sufficiently  to  maintain  its 
opacity.  The  cooling  rate  is  therefore  reduced,  the  development  is  set  back 
considerably,  and  the  cloud  does  not  manage  to  reach  the  ground  in  3  hours. 
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Figure  2.23  Case  5,  T  =  3  hours 
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Figure  2.25  Case  6,  Cold  Water,  T  =  3  hours 
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Thus,  case  6  is  not  a  simple  status  lowering  event  dependent  only  on  cloud-top 

cooling  rate,  but  is  also  strongly  affected  by  the  mixing  of  the  initial 
profile. 


Summary  and  Conclusions 

The  conclusions  which  we  are  able  to  draw  from  these  model  comparisons 
are  limited  because  we  only  have  temperature  and  humidity  profiles  at  a  few 
isolated  stations.  Discrepancies  between  model  predictions  and  verification 
data  may  be  due  either  to  modeling  deficiencies  or  to  the  neglect  of  all 
horizontal  advective  effects.  We  believe  it  is  possible  to  make  some 
conclusions  about  the  radiation  schemes. 

In  general,  the  long  wave  cooling  seems  to  be  reasonably  well  represented 
by  our  convolution  integral  and  in  most  cases  even  by  the  approximate 
convolution.  Case  4  is  the  most  straightforward  nocturnal  cooling  case,  and 
this  is  the  most  accurately  predicted  case  using  the  new  finite-differencing 
scheme  for  the  radiative  fluxes.  Cases  1,  3,  and  6  all  involve  shortwave 
solar  radiation  for  at  least  part  of  the  time,  so  that  interpretation  is  more 
difficult.  It  appeared  that  the  shortwave  heating  was  too  strong  in  Case  6, 
where  the  solar  angle  was  low.  However,  Case  1  is  similar  in  this  respect, 
and  also  produces  large  shortwave  heating,  but  still  remains  colder  than  the 
observations.  Case  3  shows  good  agreement  with  the  temperature  rise  during 
the  day,  but  this  may  be  inconclusive  since  the  cloud  continues  to  absorb 
strongly  until  it  evaporates;  thus,  if  the  absorption  was  smaller,  the  cloud 
would  presumably  remain  longer  but  could  still  absorb  the  same  amount  of  heat. 
The  heating  rates  in  clear  air  are  unaffected  by  the  scattering  formulation, 
therefore  it  is  possible  that  the  predictions  could  be  relatively  insensitive 
to  the  details  of  the  radiation  scheme.  Comparisons  with  data  cases  which 
include  direct  measurement  of  radiative  fluxes  are  necessary  if  we  are  to 
determine  the  reliability  of  the  radiation  schemes. 

We  cannot  really  conclude  much  about  the  turbulence  modeling  on  the  basis 
of  these  runs.  Although  it  appeared  that  Case  6  was  inaccurately  predicted 
partly  due  to  rapid  mixing  of  the  initial  profile,  we  do  not  know  whether  this 
was  a  model  problem  or  whether  neglected  factors  such  as  horizontal  advection 
were  important.  This  remark  also  applies  to  Case  5,  where  the  transport  of 
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humidity  was  not  predicted  during  the  first  hour,  yet  the  transport  of 
humidity  and  temperature  are  so  dissimilar  in  this  case  that  the  presence  of 
horizontal  advection  seems  very  likely. 

Regarding  the  modifications  made  to  the  model,  we  can  say  something  about 
their  relative  importance.  The  inclusion  of  the  full  atmospheric  temperature 
profiles  did  not  make  a  profound  difference  to  the  results.  Our  original 
estimates  of  the  long  wave  flux  from  aloft  were  within  20-30  Wm”^  in  most  of 
the  cases,  so  that  this  change  resulted  in  a  1°C  temperature  change  at  most. 

The  removal  of  long  wave  reflection  at  the  surface  had  a  large  effect  on 
the  stable  cases,  2  and  5.  Although  we  still  do  not  predict  fog  formation, 
the  spurious  heating  is  absent,  and  the  temperatures  are  much  closer  to  the 
observations. 

The  change  in  the  finite-difference  scheme  only  affected  Case  4,  where 
the  entrainment  was  reduced  to  almost  zero.  It  is  strange  that  none  of  the 
other  cases  were  affected  by  this  numerical  problem,  but  we  have  no 
explanation  for  their  insensitivity. 

Finally,  the  full  convolution  integral  only  affected  Case  6 
significantly,  where  it  contributed  to  the  reduced  cooling.  This  is 
apparently  because  the  approximation  for  the  convolution  gives  a  reasonable 
representation  of  an  optically  thick  cloud,  which  most  cases  are,  but  gives 
excessive  cooling  for  thin  clouds. 
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3.  ON  THE  ONE -DIMENSIONAL  PARAMETERIZATION  OF  WAVE -TURBULENT  INTERACTIONS 

Introduction 

One  of  the  principal  mechanisms  controlling  the  development  of  the 
planetary  boundary  layer  under  well-mixed  conditions  is  the  rate  of 
entrainment  of  heat  and  momentum  from  the  free  atmosphere  above  the  turbulent 
mixing  region.  This  is  due  to  the  fact  that  the  boundary  layer  eddies  do 
produce  a  well -mixed  region,  and  therefore  it  is  possible  to  describe  many  of 
the  gross  features  such  as  boundary-layer  depth,  mean  temperature,  and  mean 
velocity  by  use  of  a  'slab'  model  which  only  requires  knowledge  of  fluxes  at 
the  top  and  bottom  of  the  slab.  The  'slab'  model  does  not  predict  details  of 
the  turbulence  variations  across  the  layer,  of  course;  but  the  gross  dynamics 
of  the  actual  physical  system,  and  therefore  any  sophisticated  mathematical 
model  of  it,  are  also  controlled  by  these  surface  and  entrainment  fluxes. 

There  are  several  mechanisms  which  are  responsible  for  fluxes  at  the  top 
of  the  boundary  layer.  The  relative  importance  of  the  various  mechanisms  will 
be  determined  both  by  the  details  of  the  profiles  in  the  atmosphere  above  the 
boundary  layer,  and  by  the  boundary-layer  turbulence  itself  which  is  the 
driving  force  for  the  entrainment.  The  atmosphere  above  the  boundary  layer  is 
generally  stably-stratified  and  thus  inhibits  vertical  mixing,  since  the 
turbulence  has  to  expend  its  kinetic  energy  in  order  to  provide  the  increase 
in  potential  energy  required  to  transport  the  overlying  warm  air  down  into  the 
mixed  region. 

Since  the  air  from  just  above  the  boundary  layer  is  rapidly  cooled  to  the 

average  mixed-layer  temperature,  there  is  usually  a  significant  temperature 

change  across  a  relatively  short  distance  at  the  top  of  the  layer;  this  is 

the  so-called  "capping  inversion".  The  first  type  of  entrainment  mechanism  is 

due  to  rapidly-rising  thermals  which,  under  convective  conditions,  have 

sufficient  vertical  momentum  upon  reaching  the  inversion  to  continue  rising 

for  some  distance.  The  thermal  then  spreads  out  in  the  horizontal  and  falls 

back  into  the  boundary  layer,  entraining  some  of  the  free  atmosphere  as  it 

does  so.  This  mechanism  is  largely  controlled  by  the  vertical  velocities  in 

the  thermals,  which  in  turn  depend  strongly  on  the  surface  heat  flux.  A 
.  ,  .  .  20  21 
simple  parameterization,  e.g.,  Carson  ,  Stull  ,  sets  the  entrained  heat  flux 
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proportional  to  the  surface  heat  flux  under  conditions  of  strong  convection. 


A  second  entrainment  mechanism  which  will  be  more  important  under  more 
neutral  conditions,  i.e.,  higher  wind  speeds  or  smaller  surface  heat  fluxes, 
is  stratified  shear  instability  of  the  inversion  itself.  Since  the  inversion 
is  continually  perturbed  by  turbulent  eddies  from  below  and  possibly  by  wave 
disturbances  from  above,  we  have  a  situation  where  the  inversion  can  break 
down  locally  into  relatively  small  scale  turbulent  patches  which  entrain  the 
free  atmospheric  air. 

We  are  interested  here  in  this  latter  mechanism,  which  seems  more 
difficult  to  parameterize,  in  view  of  the  fact  that  it  can  be  affected  by  both 
local  conditions  at  the  inversion  and  the  perturbations  produced  by  the 
turbulence  in  the  boundary  layer.  In  attempting  to  gain  some  insight  into  the 
physical  processes,  we  have  first  studied  idealized  stratified  shear 
instabilities  and  their  development  through  turbulent  mixing  to  a  final  mixed 
state.  The  response  of  the  idealized  inversion  layers  to  imposed  disturbances 
has  provided  a  basis  for  understanding  the  more  general  atmospheric  problem. 
In  the  next  section,  we  present  a  summary  of  the  numerical  experiments 
performed  with  the  two-dimensional,  second-order  closure  model  on  simple 
stratified  shear  profiles.  Finally,  in  Sections  3.3  and  3.4,  we  indicate  the 
application  of  these  studies  in  the  development  of  a  simpler  model  of  the 
entrainment  process  as  it  occurs  at  the  capping  inversion  of  the  atmospheric 
boundary  layer. 


Summary  of  Numerical  Results 

2 

In  a  previous  report  (Lewellen,  et  al .  )  the  application  of  the 
two-dimensional,  A.R.A.P.  second-order  closure  model  to  the  problem  of  Kelvin- 
Helmholtz  billow  growth  and  breakdown  into  turbulence  was  described  in  detail. 
It  was  demonstrated  that  the  closure  model  has  considerable  merit  as  a  means 
of  investigating  the  detailed  dynamics  of  the  turbulent  breakdown  process. 
The  results  have  been  extended  to  investigate  the  sensitivity  of  the  mixing 
process  and  final  state  to  changes  in  the  initial  profiles  and  in  the  imposed 
initial  disturbance.  We  may  say  here  that  it  appears  there  is  one  very  simple 
result  -  namely,  that  the  final  state  after  turbulent  mixing  is  the  same  for  a 
wide  range  of  initial  conditions.  This  mixed  layer  has  a  nearly  constant 
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Richardson  number  roughly  between  0.35  and  0.40.  We  now  proceed  to  discuss 
these  numerical  results  in  more  detail. 


Previous  Results 

2 

The  results  of  Lewellen,  et  al .  indicated  that  the  final  Richardson 
number  in  the  mixed  layer  varied  from  about  0.3  to  0.4  as  the  initial 
Richardson  number  varied  from  0.1  to  0.2.  Sensitivity  tests  showed  that  the 
most  sensitive  initial  parameter  was  the  turbulence  length  scale,  which  needed 
to  be  set  at  some  large  fraction  of  the  initial  shear  layer  thickness.  The 
initial  profile  for  these  studies  was  taken  to  be 


u  =Au  tanhz/<S  ,  T  =AT  tanh  z/6 


where  u  is  the  horizontal  velocity,  and  T  is  the  temperature  perturbation 
(which  is  proportional  to  density  perturbation  in  the  Boussinesq  equations 
which  are  used  in  the  model).  The  numerical  results  proved  relatively 
insensitive  to  variations  in  initial  vorticity  perturbation  amplitude  or 
initial  turbulence  energy  level. 

The  main  gross  dynamical  quantities  presented  were  the  total  large-scale 
roll  energy 


(u^  +  w^ )  dxdz. 


(3.2) 


and  the  small  scale  turbulence  energy 


(3.3) 
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where 


q2  =  u'2  +  v'2  +  w'2  .  (3.4) 

The  development  of  these  dimensionless  energies  for  the  cases  with 
Ri  =  0.1  and  0.2  are  shown  in  Figure  3.1.  Here  Ri  =  ( g/TQ)AT6 /A ,  where  g  is 
the  gravitational  acceleration,  and  Tq  is  the  mean  temperature  of  the  layer, 
i.e.  Ri  is  the  minimum  Richardson  number  in  the  initial  profile.  The 
dimensionless  time  r  is  defined  as  t  =  t  (gAT/ToAu).  Figure  3.1  shows  the 
initial  growth  of  the  large-scale  roll  instability  and  its  subsequent 
breakdown  and  generation  of  small  scale  turbulence  energy.  At  the  end  of  the 
integration,  the  large-scale  perturbation  has  completely  decayed  so  that  we 
have  a  horizontally  homogeneous  flow  again,  and  the  turbulence  is  also 
decayi ng. 

The  initial  and  final  Richardson  numbers  for  those  two  runs  are  shown  in 
Figure  3.2.  This  figure  shows  the  mixing  effect  of  the  instability.  In  both 
cases,  the  Richardson  number  of  the  mixed  layer  is  almost  constant,  although 
the  actual  value  is  somewhat  higher  for  the  higher  initial  Ri.  For  Ri  =  0.2, 
the  final  value  is  close  to  0.4. 

Variations  in  Initial  Profile  Shape 

In  a  real  billow  turbulence  event,  the  initial  profiles  will  not  be 
precisely  identical  hyperbolic  tangents,  and  we  therefore  need  some 
information  about  the  dependence  of  the  phenomenon  on  the  initial  profile 
shapes.  One  special  feature  of  the  previously  used  profiles  was  that  the 
temperature  and  velocity  both  had  the  same  vertical  length  scale.  Hazel^^ 
performed  linear,  inviscid  stability  analyses  of  a  wide  range  of  profile 
shapes,  and  his  results  seem  to  show  that  the  most  profound  changes  in 
stability  characteristics  are  caused  by  making  the  temperature  profile  change 
across  a  thinner  layer  than  the  velocity,  i.e.,  u  =  Au  tanh  z/6 

T  =  AT  tanh  (3.5) 

with  R  >  1.  As  R  increases,  the  initial  Richardson  number  profile  changes 
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Figur'  3.1:  Evolution  of  large  scale  (EK)  and  small  scale  (EQ)  eddy 

kinetic  energies. 


60 


CJ 


o 

CVJ 


00 


o 


CVJ 

o 

II 


41 


Figure  3.1(b);  Ri 
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Figure  3.2;  Initial  and  final  Richardson  number  profiles. 

The  initial  profile  is  denoted  by  a  solid  line. 
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’>3ure  3.2(b):  Ri 


from  a  single  minimum  at  z  =  0  to  a  local  maximum  at  z  =  0  when  R  >  /2  with 
minima  above  and  below  the  centerline.  When  R  >  2,  z  =  0  becomes  a  global 
maximum  and  the  Richardson  number  tends  monotonical  ly  to  *  zero  as  z  . 
Detailed  analyses  of  the  particular  case  R  =  5  by  Hazel  shows  that  when  Ri*R 
(J  in  his  notation)  is  larger  than  0.25,  there  are  two  unstable  stationary 
modes,  and  when  Ri  >  0.37,  the  unstable  modes  are  no  longer  stationary.  Hazel 
also  found  that  small  changes  in  the  profile  shape  (from  hyperbolic  tangent  to 
error  function)  or  small  asymmetries  about  z  =  0  did  not  materially  alter  the 
linear  stability  characteristics.  We  have  therefore  concentrated  on  variation 
in  the  thickness  ratio,  R,  in  our  studies  here. 

Figure  3.3  shows  the  initial  and  final  Richardson  number  profiles  for  a 
case  with  R  =  1.7.  The  initial  profile  has  a  local  maximum  at  z  =  0  as 
mentioned  previously,  with  secondary  minima  at  z  =  1.4.  The  maximum  at 

z  =  0  is  0.22  and  the  minimum  is  roughly  0.15.  At  a  time  t  =  24,  the 

turbulence  has  largely  decayed,  and  the  Richardson  number  profile  shows  a 

minimum  value  of  about  0.37  and  an  average  value  of  about  0.4  across  the  mixed 
layer.  The  evolution  of  the  kinetic  energies  is  shown  in  Figure  3.4.  There 

is  a  primary  breaking  event  at  x  =  3  which  generates  most  of  the  turbulent 
energy,  but  this  is  followed  by  a  secondary  ’  break  at  x  =  6  which  gives  a 

further  increase  in  EQ.  Examination  of  the  temperature  contour  patterns  shows 
that  the  primary  billow  is  the  symmetric  mode  centered  on  z  =  0,  as  can  be 
seen  in  Figure  3.5  at  x  =  2.8  and  x  =  4.3.  This  is  presumably  because  the 
inflection  point  at  z  =  0  is  the  major  source  of  instability,  and  since  the 
Richardson  number  at  z  =  0  is  only  0.22  initially,  the  fastest  growing  model 
is  centered  there.  However,  after  this  initial  system  roll-up  and  turbulent 
breaking,  a  secondary  mode  appears.  This  mode  is  evident  at  x  =  5.8  and 

X  =  7.3,  and  seems  to  be  associated  with  the  secondary  minima  in  the  initial 

Richardson  number  profile.  The  secondary  mode  also  rolls  up  the  remaining 

vortex  sheet  above  and  below  the  primary  bil.low,  and  generates  its  own 
small  scale  turbulence.  Thus  the  billow  event  is  more  complicated  than  the 
previously  studied  cases,  but  the  final  result  of  the  mixing  process  is  to 
produce  a  layer  with  a  Richardson  number  of  about  0.4. 
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Figure  3.3:  Initial  and  final  Richardson  number  profiles  for 
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Figure  3.4:  Evolution  of  kinetic  energies  for  the  case  with  R 
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Isopleths  of  dimensionless  temperature,  (T-To)/aT, 
for  the  case  with  R  =  1.7  at 
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Figure  3.5(b); 
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Figure  3.5(c): 
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Figure  3.5(d): 


When  R  =  2.5,  the  Richardson  number  profile  has  a  maximum  at  z  =  0  as  can 
be  seen  in  the  initial  profile  in  Figure  3.6.  Also  shown  in  this  figure  are 
the  late  time  profiles  from  two  numerical  integrations  using  the  given  initial 
profile.  In  the  first  run,  a  short  domain  of  length  126  was  used  with  an 
initial  perturbation  of  vorticity  amplitude.  In  the  second,  a  domain  of 
length  246  and  an  initial  perturbation  of  the  isotherms  was  applied.  The 
evolution  of  the  energies  from  the  runs  is  shown  in  Figure  3.7.  It  is  obvious 
that  very  different  modes  are  excited  in  the  two  integrations.  Temperature 
contours  at  two  times  from  each  run  are  shown  in  Figures  3.8  and  3.9.  In  the 
short  domain,  only  a  very  weak  symmetric  disturbance  is  excited,  but  there  is 
enough  turbulent  growth  in  the  low  Richardson  number  regions  to  produce  the 
required  mixing  and  stabilize  the  profile.  In  the  long  domain,  the  initial 
wave  does  not  immediately  decay  in  amplitude,  but  feeds  energy  into  a 
large-scale  mode  which  does  produce  a  large  disturbance  and  convective 
instabilities  in  the  rolled-up  vortex  cores.  Turbulence  levels  are  much 
higher,  although  the  time  scales  are  also  much  longer.  However,  the  final 
result  is,  in  both  cases,  a  mixed-layer  with  a  Richardson  number  of  roughly 
0.4.  We  do  not  wish  to  dwell  on  the  details  of  these  integrations,  since 
neither  are  ideal  examples;  for  example,  the  long  domain  case  is  probably 
significantly  affected  by  the  upper  and  lower  boundary  conditions.  However, 
these  factors  do  not  seem  to  affect  the  overall  mixing  effect  of  the 
instability. 

Larger  Initial  Richardson  Numbers 

Finally,  and  briefly,  we  show  the  results  from  initial  Richardson 
numbers  closer  to  the  critical  value  of  0.25.  This  has  relevance  to  the 
problem  of  a  slowly  decreasing  Richardson  number,  as  might  be  produced  in  the 
atmosphere  by  some  large-scale  flow  feature.  In  this  case,  the  shear  layers 
would  presumably  begin  to  roll  up  before  the  Richardson  number  had  time  to 
fal 1  much  bel ow  0.25. 

Using  the  simple  hyperbolic  tangent  profiles,  and  Ri  =  0.23,  an 

integration  was  made  with  a  domain  of  length  116.  The  energy  evolution  is 
shown  in  Figure  3.10,  and  the  Richardson  number  at  t  =  8  is  shown  in 
Figure  3.11.  Although  the  mixing  is  not  completed,  it  is  clear  that  there  has 
been  the  usual  breaking  event,  and  the  final  Richardson  number  will  be  close 


51 


a 

CO 

<T3 

U 


L_ 

o 


CO 

OJ 


o 

e 

Q. 

c  c 

£- 

•r—  •r— 

O) 

^  <T3 

-O 

E 

E  E 
O  O 

3 

"O  “O 

C 

CJ)  4-> 

C 

C 

o 

o  o 

CO 

1 —  j:; 

T3 

CO 

C_ 

c 

c: 

jr 

‘r~ 

o 

f— 

•r“ 

ru  . — 

or 

C  03 

r— 

•«“  e 

rXS 

4-  T- 

^t3 

•r- 

4- 

C 

-f-> 

•r— 

4- 

C 

X5 

c 


*“  LT) 

ro  • 

•r-  OJ 
4J 

•f-  H 

c 

or 


ro 

0) 

c. 

a 

C7» 


52 


0.60 


53 


54 


4.80 


O 

CD 

Csl 


O 

CM 


o 

CD 

cj5 


L- 

o 


~o 

c 

<T3 


un 

OJ 

II 

cn 


t- 

o 


O) 

u 

o 

13 

o 

+-> 

<v 

00 

C- 

OJ 

Q_ 

e 

<v 

+j 

o 

CO 

to 

CO 

• 

OJ 

CO 

c 

o 

CO 

c 

<3J 

o 

E 

00 

•»— 

• 

•O 

4- 

O 

CO 

jr 

o 

CM 

cu  c 

• 

Q.  rc 

CO 

O  E 
CO  o 
1-^  XJ 

fX3 

> 

<v 


L- 

3 

o 

4-> 

CT 

O 

o 

to 


II 


O 

CO 


o 

o 


00 

CO 

a’ 

t- 


O 

o 

o 

o 

o 

cdo 

CM 

CO 

o 

CO 

CM 

CO 

• 

• 

• 

• 

■ 

m 

CO 

o 

•o 

M 

1 

CO 

1 

55 


4.80 


o 

CO 

CVJ 


o 

Csl 


o 

to 

o5 


o 

o 

00 


o 

to 


o 

oo 

■ 


o 

CM 

CO 


o 

to 


o 

o 


o 

O 

o 

o 

o 

oo 

CM 

to 

o 

to 

CM 

00 

CO 

o 

CO 

1 

J 

I 

tM 


K3 


LO 


It 

h- 


56 


Figure  3.8(b); 
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Figure  3.9;  As  Figure  3.8  but  with  the  long  domain. 
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Fic/ire  3.9(c): 
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Evolution  of  kinetic  energies  for  the  case  with  Ri  =  0.23 
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to  0.4.  In  fact,  the  temperature  contour  plots  show  a  breaking  billow  very 
similar  to  the  Ri  =  0.2  case. 

This  result  raises  the  question  of  whether  the  instability  is  quenched 
when  Ri  >  0.25.  Thus,  a  case  with  Ri  =  0.27  was  run  with  all  other  parameters 
identical  to  the  previous  run.  The  energy  evolution  and  final  Richardson 
number  profiles  are  shown  in  Figures  3.12  and  3.13.  Once  again,  there  is  a 
vortex  roll-up  and  collapse  giving  a  final  Richardson  number  of  about  0.4. 
Clearly,  there  is  a  finite  amplitude  instability  of  the  shear  layer  profile 
for  Ri  =  0.27  (linear  theory  predicts  stability),  and  our  initial  perturbation 

of  amplitude  10%  of  the  background  vorticity  is  sufficiently  large  to  trigger'" 

it. 


We  have  not  performed  extensive  studies  of  the  dependence  of  the  flow  on 
initial  perturbations  or  length  of  integration  domain,  and  so  we  can  only  say 
that  the  triggering  of  atmospheric  billow  events  is  likely  to  occur  at 
Ri  =0.25  and  will  depend  on  a  finite  amplitude  perturbation,  unless  the 
Richardson  number  is  reduced  very  quickly  into  the  linear  instability  regime. 

One-Dimensional  Calculations  of  the  Kelvin-Helmholtz  Instability 

Our  two-dimensional  model  of  the  stratified  shear  instability  has 
demonstrated  the  mechanism  for  billow  breakdown  into  small  scale  turbulence, 
and  also  the  ability  of  the  second-order  closure  scheme  to  model  this  process. 
In  these  detailed  integrations,  attention  was  focussed  on  the  individual 
billow  event,  and  the  initial  coherent  vortex  roll-up  stage  was  calculated 
explicitly.  In  a  larger  scale  problem  where  individual  billows  are  no  longer 
resolved,  the  closure  model  has  to  describe  the  vortex  roll-up  stage  as  well 
as  the  vortex  breakdown.  In  this  section,  we  discuss  the  problems  involved  in 
such  a  calculation. 

The  first  point  to  emphasize  is  that  the  initial  stage  of  the  instability 
is  linear,  giving  exponential  growth  of  a  coherent  wave.  We  immediately 
encountered  problems  with  the  A.R.A.P.  model  because  it  does  not  predict  any 
linear  exponential  growth  for  Richardson  number  greater  than  or  equal  to  zero. 
Exponential  growth  is  only  achieved  when  the  nonlinear  terms,  i.e.,  the 
return-to-i sotropy  term,  is  significant,  and  the  growth  rate  in  this  regime  is 
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almost  an  order  of  magnitude  slower  than  the  linear  stability  analysis 
prediction  for  the  hyperbolic  tangent  profiles. 

The  deficiency  of  the  A.R.A.P.  model  is  easily  traced  to  the  absence  of 
"rapid"  pressure  terms.  If  we  assume  that  the  turbulence  energy,  q  ,  is  very 
small,  so  that  q/A  «  dU/dz  or  N,  the  Brunt  Vaisaila  frequency,  then  only  the 
production  terms  are  important.  For  the  shear  layer,  the  turbulence  needs  to 
develop  a  considerable  degree  of  anistopy  through  shear  production  before  the 
Rotta  term  can  provide  the  positive  feedback  which  we  require  for  exponential 
growth.  We  should  point  out  here  that  this  discussion  does  not  invalidate  the 
two-dimensional  integrations,  where  we  required  the  small  scale  turbulence  in 
the  billow  to  grow  exponentially  from  small  levels  to  model  the 
three-dimensional  secondary  instability.  The  mechanism  for  the  latter  was 
shown  to  be  the  convective  instability  produced  by  overturning  the  temperature 
gradient;  the  A.R.A.P.  model  does  predict  exponential  growth  of  this  mode, 
because  there  is  a  direct  feedback  from  the  heat  flux,  into  the  vertical 

energy  component,  w^. 

The  remedy  in  the  case  of  the  shear  layer  appears  to  be  the  inclusion  of 

23 

"rapid"  pressure  terms,  as  advocated  by  Hanjalic  and  Launder  ,  and  Lumley  and 

.24 

Khajeh-Nouri  .  These  terms  redistribute  the  production  of  energy  between  the 
tensor  components  instantaneously,  and  arise  from  the  pressure  fluctuations 
driven  by  the  mean  velocity  gradient  and  the  buoyancy  fluctuations.  There  is, 
however,  still  a  problem  in  that  the  model  terms  suggested  by  e.g.,  Gibson  and 
Launder  do  not  give  a  critical  Richardson  number  of  0.25. 

Briefly,  Gibson  and  Launder  write  the  Reynolds-averaged  equation  as 
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where 
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J  ax 


Gij  =  gi  Uj0  +  gj  Ui0  .  G^-  =  g^  ©2  (3.7) 

where  (j),  D,  and  e  represent  pressure  correlations,  triple  correlations  and 
dissipation,  respectively.  The  A.R.A.P.  model  sets 


♦u  “  ♦i-j*  -  -  A  ^  ;  *1  "  ♦S'*  =  -  A  a  u,e  (3.8) 


2  -  (1  ) 

where  =  u.jU.j,  and  A  is  the  turbulence  length-scale.  ')>jj  is  the  Rotta 

term. 


25 


Gibson  and  Launder  recommend  the  inclusion  of  extra  pressure  terms: 


and 


J2) 

<j)i  =  -  CZq  Pi 


(3) 


C30  Gi 


(3.9) 


(3.10) 


where  0^=0. 6,  Cg=0.5,  C2^=C2^=0.33.  Thep^coefficients  also  involve  a  change 
in  the  coefficient  on  the  Rotta  term,  from  unity  to  0.45.  Note  that  the 
production  term  -u-jUj  30/axj  does  not  give  rise  to  any  rapid  pressure  term, 
because  this  term  does  not  come  from  multiplying  the  momentum  equations,  as 
discussed  by  Launder  . 
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We  have  used  the  above  model  for  the  rapid  terms  (but  used  0.375  for  the 
Rotta  coefficient  in  line  with  Launder,  Reece,  and  Rodi^^,  since  this  gave 
better  agreement  in  the  neutral  surface  layer)  in  the  atmospheric  surface 
layer,  and  obtained  results  comparable  with  the  A.R.A.P.  model.  (Note  that 
these  results  were  obtained  without  the  complicated  wall  functions  of  Gibson 
and  Launder.)  Thus,  the  rapid  terms  do  not  degrade  the  model  performance  in 
this  regime;  ^  although  it  must  be  admitted  that  the  inclusion  of  four  extra 
empirical  constants  does  not  significantly  improve  the  results,  either. 
However,  we  do  have  rapid  terms  which  can  promote  exponential  growth  in  the 
shear  layer. 

Unfortunately,  the  critical  Richardson  number  for  this  model  in  the 
linear  regime  is  0.15,  which  is  significantly  smaller  than  the  actual  value  of 
1/4.  Adjusting  the  constants  to  move  the  critical  value  to  1/4  degraded  the 
surface  layer  performance  of  the  model  quite  seriously.  The  remedy  we  found 

was^o  include  a  rapid  term  proportional  to  the  scalar  flux  production  term, 

-u^Uj  30/3Xj.  If  we  apply  a  coefficient  of  1/3  to  this  term,  i.e.,  the  entire 

production  of  u^e  is  reduced  by  1/3,  and  set  c^=c^=0.6,  then  this  simple  model 

has  a  critical  Richardson  number  of  0.25,  and  produces  very  similar  results  to 
the  A.R.A.P.  model  in  the  surface  layer.  The  linear  growth  rates  predicted  by 
the  model  are  actually  very  close  to  the  hyperbolic  tangent  profile  growth 
rates  calculated  by  Hazel  ,  as  can  be  seen  from  Figure  3.14. 

A  calculation  for  the  tanh  profile  shows  reasonable  agreement  with  the 
two-dimensional  calculations,  provided  we  fix  the  length  at  a  reasonable 
fraction  of  the  shear  layer  thickness.  Figure  3.15  shows  the  evolution  of  the 
Integrated  turbulence  energy  for  a  case  with  Ri  =  0.1  initially.  The  time 
scales  and  maximum  values  are  close  to  the  two-dimensional  values,  as  is  the 
final  Richardson  number  as  shown  in  Figure  3.16. 

The  problem  cannot  be  claimed  to  be  entirely  solved,  however,  because  in 
the  case  of  a  passive  scalar,  we  would  not  want  to  include  any  rapid  term 

corresponding  to  30/3Xj  in  the  scalar  flux  equation.  Analytical 

solutions  for  the  initial  rate  of  diffusion  from  a  point  source  show  that 
there  is  no  such  rapid  term.  We  are  therefore  unable  to  present  a  truly 
invariant  model  which  adequately  describes  both  the  initial  growth  of 

turbulence  in  a  shear  layer  and  the  initial  rate  of  diffusion  from  a  scalar 
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"  Model  Result 

O — —O  Hyperbolic  Tangent  Profile  (analytic) 


Figure  3.14:  Dimensionless  growth  rates  for  stratified  shear  layer. 

Perturbations  are  proportional  to  exp(at),  and  S  is  the 
velocity  shear. 

_ closure  model,  homogeneous  shear 

-  linear  analysis  for  hyperooiic 

tangent  profile  (Hazel,  1972) 
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Figure  3.15:  Evolution  of  total  energy,  EQ,  for  one-dimensi oriu  i  closure 

model  integration.  Horizontal  length  of  the  dom,r  vas  taken 
to  be  165  for  direct  comparison  v/ith  two-diinensionu  •  ..ases. 
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closure  model  integration. 


point  source. 


Development  of  Parameterizations 


The  knowledge  of  the  detailed  dynamics  of  the  Kel vin-Helmhol tz  breaking 
process  can  be  used  to  develop  both  simple  parameterizations  of  the  billow- 
induced  mixing  itself,  and  also  to  study  more  complex  entrainment  processes 

when  the  small  scale  shear  instability  is  a  secondary  mechanism. 

In  the  former  case,  we  are  concerned  with  Kel vin-Helmholtz  billows  as  a 

dominant  feature  of  the  inversion.  From  our  detailed  studies,  we  may  assert 

that  whenever  the  inversion  rolls  up  and  breaks  down  into  turbulence,  the 

final  mixed  state  will  be  a  shear  layer  with  a  Richardson  number  of  about  0.4. 

Thus,  if  we  can  estimate  the  initial  state  before  the  roll-up  we  can  calculate 

the  total  entrainment;  if  we  can  further  estimate  the  frequency  of  the 

breaking  events,  then  we  have  a  value  for  the  average  entrainment  rate  at  the 

2  8 

interface.  In  Stull's  parameterization  of  shear-induced  mixing  at  the 
inversion,  the  entrainment  velocity  is  expressed  in  terms  of  inversion 

parameters  only,  i.e.,  thickness,  6,  velocity  jump,  aU,  and  temperature  jump, 
aT.  This  may  be  an  oversimplification,  since  the  breaking  events  must  be 

triggered  by  external  perturbation  or  there  will  be  one  event  which  produces  a 
Richardson  number  of  0.4  and  no  further  activity.  Thus,  one  aspect  of  the 
problem  which  deserves  further  attention  is  the  form  of  the  perturbations  at 
the  interface  and  their  dependence  on  the  boundary  layer  turbulence.  For  the 
moment,  we  can  give  an  illustration  of  the  parameterization  method  by  making 
assumptions  about  the  interface  perturbations.  The  temperature  equation 
states  that 


where  T  is  the  temperature,  and  H  is  the  heat  flux.  Thus,  if  we  define 


(3.12) 
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where  the  integral  is  carried  out  over  the  entire  billow  event,  then 

1^  =  To  (z)  -  Tf(z).  (3.13) 

where  Tq  and  T|:  are  the  initial  and  final  temperature  profiles. 

For  the  sake  of  simplicity,  let  us  assume  piecewise  linear  profiles  so 

that 


r 


AT/2 


^0,f  ~  ^  ^^^/^0,f 

-AT/2 


^  <  ^o,f/2 


z  <  -  Jlo,f/2 


(3.14) 


where  Aq*  ^f  initial  and  final  shear  layer  thickness.  If  we  assume 

that  the  velocity  profile  has  the  same  shape,  then  the  initial  Richardson 
number.  Rig  =  (g/To)^T^-o/^u^,  and  the  final  value. 


ATJlf 

Au^ 


(3.15) 


Now  we  know  Ri^  =  0.4.  Thus 


(3.16) 


For  definiteness,  let  us  assume  that  Rig  =  0.133,  so  that  if  =  3jIo;  we  can 
now  calculate  the  integrated  heat  flux.  To  obtain  the  maximum  value,  F(0),  we 
have 


F(0).-/  {Tj,(z)  -  Tf{z))dz 

0 
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(3.17) 


To  obtain  the  average  heat  flux,  H,  at  the  inversion,  we  now  need  to 
estimate  the  frequency  of  those  breaking  events;  then 


H  =  =  ^0^"^ 


(3.18) 


where  t  is  the  period  between  events. 

As  we  have  already  stated,  t  will  depend  on  the  boundary  layer 
turbulence,  but  a  lower  limit  on  t  is  given  by  the  time  scale  of  the  breaking 
event.  This  is  because  the  layer  must  recover  from  the  turbulent  breakdown, 
otherwise  we  do  not  have  distinct  events.  From  our  detailed  calculations,  the 
time  scale  is  on  the  order  of  lOa/Au,  where  a  is  the  wavelength  of  the 
disturbance,  which  is  roughly  10£q.  Thus  a  minimum  value  for  t  is  100£q/Au. 
Therefore 


£-  AT  Au 


(3.19) 


To  obtain  our  entrainment  velocity,  Wg,  we  divide  H  by  AT  giving  Wg  <  2.5  x 
10"3  Au. 
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Stull  parameterizes  the  entrainment  velocity  as 


(3.20) 
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where  Rig  is  the  bulk  Richardson  number  of  the  inversion,  and  Ag  is  an 
empirical  constant  which  is  derived  from  experimental  results  to  be  10"^.  If 
we  take  Rig  =  0.4  as  a  typical  value  of  the  inversion,  we  obtain  precisely  the 

same  estimate  as  our  crude  upper  bound.  The  exact  agreement  is  obviously 

fortuitous,  since  we  only  made  rough  estimates  of  the  effects  of  the  boundary 
layer  turbulence  -  but  is  nevertheless  encouraging.  The  experiments  used  in 
deriving  the  value  for  A  were  explicitly  concerned  with  the  shear-induced 
instability,  so  we  would  expect  the  value  to  be  near  our  upper  bound.  More 
work  is  required  to  determine  the  dependence  of  A^  on  the  boundary  layer 

turbul ence. 

The  second  application  of  the  detailed  study  is  in  the  parameterization 
of  shear  instability  as  a  secondary  mechanism;  for  example,  in  the  case  of 
entrainment  by  penetrative  convection.  These  processes  may  be  studied  using  a 
numerical  model  of  the  large  scale  features  in  conjunction  with  a  simple  model 
of  the  small  scale  processes.  In  this  regard,  the  one-dimensional 

integrations  test  the  ability  of  the  second-order  closure  model  to  describe 
the  overall  mixing  event  without  resolving  the  detailed  billows  themselves. 
Our  integrations  show  that  this  can  be  achieved  with  reasonable  success 
provided  that  the  length  scales  can  be  adequately  described.  This  may  require 
more  effort  in  the  specification  of  the  length  scale,  for  example,  resetting 
the  length  scale  to  some  specified  fraction  of  the  shear  layer  thickness 
whenever  the  Richardson  number  falls  below  some  critical  value  and  the 
turbulence  energy  is  small.  It  does  seem  likely  that  these  small  scale 
processes  can  be  described  as  sub-grid  mixing  using  the  second-order  closure 
model,  which  will  permit  the  further  study  of  larger  scale  inversion  dynamics. 

As  noted  earlier,  in  considering  the  problem  of  entrainment  at  the  top  of 
the  atmospheric  boundary  layer,  the  small  scale  shear  instabilities  are 
generally  a  secondary  mechanism.  Having  achieved  some  understanding  of  these 
small  scale  processes,  the  next  problem  is  to  investigate  how  they  are 
triggered  by  the  turbulent  boundary  layer  eddies.  This  requires  consideration 
of  the  whole  mixing  layer,  and  there  are  two  feasible  numerical  approaches. 
One  can  choose  to  model  the  boundary  layer  eddies  in  either  two  or  three 
dimensions.  Three  dimensions  is  clearly  more  realistic  but  is  much  more 
expensive  and  demands  a  very  simple  sub-grid  closure.  The  simple  closures 
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have  proven  very  effective  in  neutral  and  convectively  unstable  flows,  but 
have  not  been  tested  sufficiently  in  stratified  cases.  On  the  other  hand, 
two-dimensional  calculations  have  reproduced  many  of  the  observed  features  of 
boundary-layer  eddies,  and  also  have  some  observational  justification  insofar 
as  longitudinal  rolls  are  quite  common  under  certain  circumstances.  The 
two-dimensional  calculations  would  also  permit  a  more  sophisticated  sub-grid 
closure,  which  would  give  more  confidence  in  the  stratified  shear  layer  at  the 
inversion. 
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4.  COMMENTS  ON  THE  EXTENSION  OF  OUR  INTEGRAL  BOUNDARY  LAYER  MODEL 
FROM  ONE-DIMENSIONAL  TO  THREE-DIMENSIONAL 

As  detailed  in  past  reports  (References  2  and  17)  we  have  spent  some  time 
attempting  to  reduce  the  grid  resolution  required  to  adequately  represent  the 
atmospheric  boundary  layer  by  imposing  integral  constraints  on  the  finite 
difference  algorithms.  The  goal  has  been  to  develop  a  hybrid 

integral -di fferential  method  which  can  be  used  with  any  number  of  vertical 
layers  from  1  to  10.  The  single  layer  representation  is  purely  integral  but 
should  provide  (at  best)  a  rough  approximation  to  the  boundary-layer  dynamics. 
The  extent  to  which  we  have  been  able  to  do  this  for  the  homogeneous  boundary 
layer  in  some  different  dynamical  situations  is  illustrated  in  Figures  3.2  to 
3.5  of  Reference  2.  Some  effort  has  been  expended  on  continuing  this  approach 
during  the  current  contract  period.  These  efforts  have  convinced  us  that  the 
most  promising  way  to  proceed  with  this  approach  is  to  combine  it  with  the 
general  problem  of  sub-grid  flux  parameterization.  This  may  be  elucidated  by 
discussing  the  problems  and  promises  of  using  this  integral  approach  for  a 
fully  three-dimensional  boundary  layer  model. 

If  the  boundary  layer  can  be  adequately  represented  by  only  a  few 

vertical  layers,  then  the  three-dimensional,  boundary-layer  problem  in  x,y,z, 

and  t  can  be  reduced  to  a  quasi-two-dimensional  problem  in  x,y,  and  t.  This 

very  attractive  possibility  has  supplied  most  of  the  impetus  for  our  integral 

modeling  efforts.  Single  layer  representations  of  the  planetary  boundary 

layer  (PBL)  for  such  a  reduced  three-dimensional  problem  have  met  with  mixed 

29 

success  in  the  literature.  The  simple  mixed  layer  models  of  Lavoie  and  by 
30 

Keyser  and  Anthes  appear  to  provide  considerable  realism  for  the  relatively 
small  computing  requirements  they  have  in  comparison  to  multi-level,  fully 

31  3233 

three-dimensional  models  such  as  Warner  et  al .  .  However,  Anthes  et  al .  ’ 

have  shown  that,  at  least,  their  version  of  a  mixed  layer  PBL  is  not  able  to 
adequately  represent  the  flow  within  the  PBL  when  horizontal  inhomogeneities 
associated  with  differential  heating  over  complex  terrain  or  across  land-water 
boundaries  is  a  dominant  mechanism. 

Anthes,  et  al.,  attribute  this  deficiency  of  the  mixed  layer  model  to 
difficulties  in  representing  the  horizontal  pressure  gradient  at  the  top  of 
their  mixed  layer.  However,  in  more  general  terms,  it  appears  to  be  a 
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reflection  of  the  difficulty  of  providing  a  simple  parameterization  which  can 
represent  a  wide  variety  of  profiles  of  the  mean  variables  within  the  PBL. 
The  sea  breeze  example,  which  they  consider,  represents  a  case  where  the 
boundary  layer  dynamics  impose  a  circulatory  flow  pattern  in  the  vicinity  of 
the  shoreline.  They  identified  the  on-shore  flow  as  occurring  within  their 
mixed  layer  PBL  and  forced  the  return  flow  to  occur  above  it.  Their  analysis 
then  shows  that  the  parameterization  of  the  return  flow  layer  is  as  important 
as  the  parameterization  of  the  mixed  layer  itself.  It  is  perhaps  more 
appropriate  to  think  of  the  return  flow  as  the  outer  part  of  the  total  PBL. 
When  viewed  in  this  way,  the  sea  breeze  involves  flow  with  a  reversal  in 
direction  of  the  velocity  within  the  PBL.  The  representation  of  reversed  flow 

profiles  necessarily  complicates  any  parameterization  of  a  single  layer  model. 

34 

The  history  of  integral  boundary  layer  models  in  aeronautics  (Schlichting  ) 
have  shown  them  to  have  quite  limited  success  under  separated  flow  conditions 
which  yield  a  flow  reversal  in  the  boundary  layer.  This  leads  us  to  be  rather 
pessimistic  about  the  prospective  adequacy  of  any  single  layer  representation 
of  the  PBL  for  general  three-dimensional  meteorological  problems.  This  still 
leaves  the  possibility  that  a  model  with  a  few  layers  can  be  quite  successful. 

The  extension  of  our  hybrid  integral  model  to  fully  three-dimensional 

situations  requires  extending  the  complete  second-order  closure  model  to 

three-dimensions  first.  Problems  associated  with  programming  algorithms  for 

computing  the  mean  flow  variables  in  existing  three-dimensional,  primitive 
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Equation  models  such  as  Pielke's  must  be  tackled  along  with  those  required 
for  computing  the  second-order  fluxes.  A  task  that  should  be  simpler,  and 
thus  of  more  immediate  utility,  is  to  use  our  turbulent  transport  model  to 
construct  algorithms  to  parameterize  the  role  of  sub-grid  turbulent  fluxes  in 
existing  mesoscale  meteorological  models.  Integral  constraints  generated  by 
integrating  the  second-order  flux  equations  over  the  resolved  grid  lengths  can 
be  very  useful  in  this  sub-grid  parameterization  role  without  the  need  to 
construct  a  completely  new  three-dimensional  model. 

As  a  simple  example  to  illustrate  the  latter  idea,  consider  the  sub-grid 
flux  parameterization  of  the  lowest  layer  of  the  flow  where  the  vertical  grid 
spacing  is  such  as  to  completely  contain  the  boundary  layer  within  the  lowest 
level.  The  flux  term  in  the  u  momentum  equation  we  wish  to  parameterize  is 
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(4.1) 


9u'w' 

3Z 


The  level  1  Equation  for  u  then  involves  some  representation  of 


(u'w'  -  u'Wq)/2az 


(4.2) 


If  the  eddy  parameterization  is  used,  then 


u'w'  =  -  e  8U/9Z 


(4.3) 


and  a  very  poor  representation  of  the  surface  shear  stress  is 
a  compensating  slip  velocity  is  permitted  at  the  surface, 
form  for  u'w'  can  be  obtained  for  the  superequilibrium 
turbulent  correlation  equations 


obtained,  unless 
A  simple  useful 
balance  of  the 


u  w  = 


/2  .2 


ly. 
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az 


(4.4) 


This  mixing  length  form  of  the  parameterization  is  quite  appropriate  in  the 
surface  layer  with  A  a  z,  but  some  integral  form  of  it  must  be  used  to  relate 
the  surface  shear  stress  to  the  velocity  at  the  first  grid  point.  If  the 
first  grid  point  is  placed  within  the  surface  layer  then 


u'w')o 


k 


(An  z^/Zg)' 


(4.5) 


where  k=0.4  and  Zq  is  the  effective  aerodynamic  roughness  of  the  surface.  A 
combination  of  Equations  4.4  and  4.5  is  quite  effective  in  Equation  4.2,  when 
the  boundary  layer  is  fairly  wel 1 -resol ved  and  there  is  no  influence  of 
stratification  or  pressure  gradient  on  the  turbulence.  However,  even  in  this 
simple  case  it  is  clear  that  Equation  4.5  needs  some  adjustment  if  the 
resolution  is  reduced  so  that  the  first  grid  point  lies  outside  the  constant 
flux  surface  layer. 
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In  sunmary,  we  believe  it  is  appropriate  to  change  the  present  emphasis 
of  our  integral  modeling  effort  from  that  aimed  at  developing  a  fully 
three-dimensional  version  of  our  second-order  closure  model,  to  that  of 
providing  support  for  our  future  work  of  parameterizing  the  role  of  sub-grid 
turbulent  fluxes  for  mesoscale  meteorological  models.  This  simple  example 
above  is  meant  to  illustrate  how  we  expect  our  past  work  on  hybrid  integral 
models  of  the  PBL  to  be  quite  helpful  in  this  new  task. 
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5.  PRELIMINARY  CONSIDERATIONS  OF  CUMULUS  PARAMETERIZATION 
BASED  ON  SECOND-ORDER  CLOSURE 

Introduction 

Some  of  the  most  important  turbulent  transport  processes  in  the  lower 

troposphere  involve  cumulus  clouds.  In  global  models  this  calls  for  the 

introduction  of  a  cumulus  parameterizations  scheme  to  represent  the  turbulent 

transport  of  humidity,  heat,  and  momentum  by  these  cumulus  clouds.  Such 
,  36  37  38, 

schemes  (Arakawa  &  Schubert  ,  Lord  ,  Kuo  )  require  a  phenomenological 

description  of  the  most  important  effects  of  clouds. 

In  what  follows  we  take  up  a  preliminary  consideration  of  the  application 
of  higher-order  closure  turbulence  methodology  to  such  flows.  We  note  that 
such  a  higher-order  closure  description  should  be  capable  of  describing  a 
range  of  atmospheric  flows  from  boundary  layers  to  mesoscale  motions  in  which 
cloudiness  and  conditional  instability  are  of  importance.  In  addition,  with 
appropriate  modifications  for  the  large  density  changes  over  the  turbulent 
macroscale,  such  a  higher-order  closure  turbulence  theory  should  naturally 
describe  the  effects  of  deep  cumulus  convection  and  hence  provide  an 

alternative  approach  to  the  problem  of  cumulus  parameterization. 

The  presence  of  both  saturated  and  unsaturated  regions  within  a  buoyant 
flow  has  at  least  two  important  influences  on  the  turbulent  description  of 
such  flows.  The  first  is  the  appearance  of  new  correlations  in  the  turbulence 
equations,  which  result  from  the  turbulent  fluctuations  of  the  clouds. 
Principal  among  these  cloud  correlations  are  the  mean  cloudiness  r,  the 
cloud  flux  w'r‘,  and  the  cloud  turbulence  energy  1/2  w'w'r',  where  r  is  a 
conditional  variable,  with  the  values  r=l  in  cloud  and  r=0  in  clear  air. 

Second,  flows  that  contain  both  saturated  and  unsaturated  regions  are  very 
likely  to  be  conditionally  unstable.  Thus,  most  of  the  large  scale  field  may 
be  virtually  devoid  of  turbulence  while  local  isolated  regions  consisting  of 
cloudy  and  non-cloudy  updrafts  and  downdrafts  may  be  in  strong  convection. 

This  latter  effect  requires  the  introduction  of  the  concept  of  an  intermittent 

turbulence  field  for  the  large  scale  flow. 
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In  Part  2  we  review  the  partly  cloudy  correlations  which  must  appear  in 
the  turbulence  equations.  In  Part  3  we  consider  second-order  closure 
turbulence  modeling  of  an  intermittent  flow.  In  Part  4  we  consider  the 
application  of  this  model  and  we  make  some  data  comparisons  and  checks  with 
the  theory  of  Parts  2  and  3. 

Second-Order  Closure  Turbulence  Equations  with  Cloud  Correlations 

We  begin  with  the  instantaneous  energy  Equation  written  in  terms  of  the 

12 

virtual  potential  temperature  6y.  (Oliver,  Lewellen,  and  Williamson  ): 


DSy/Dt  =  <5  S  +  Kg  Tw 

(5.1) 

where  w  is  the 

precipitation 

vertical  velocity,  r  the  adiabatic  lapse, 
source  term.  We  may  express  <5  and  Kg  as 

and  S  the  radiant  and 

Ks  "  1  +  (3s-l)r 

(5.2a) 

'^g  ^g'' 

(5.2b) 

with 

65  =  y/u  3g  =  a-y2/y 

and  p,  y,  and  a  are  the  functions  of  the  saturation  mixing  as  given  by 
12 

01 iver  et  al . 

The  variable  r  is  the  instantaneous  cloudiness  defined  as 

r=H  (q^)  (5.3) 

where  H  is  the  Heaviside  function  and  q^  is  the  difference  mixing  ratio 

qx  =  q-qs 

with  q^  the  saturation  mixing  ratio.  The  instantaneous  liquid  water  mixing 
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ratio  is  given  by 


qji  = 


(5.5) 


We  note  that  because  of  the  sensitive  dependence  of  radiative  transport 
on  the  presence  of  liquid  water  drops,  the  radiant  flux  divergence  S  will 
depend  upon  the  cloudiness  present  (as  well  as  the  water  vapor  in  both  clear 
and  cloudy  air).  We  do  not  take  up  the  procedure  of  treating  the  ensemble 
averaging  of  S  here.  We  merely  indicate  its  presence  through  the  source  term 
Z  defined  as 


=  [l+(3s-l)r]S 


(5.6) 


We  may  express  Equation  (5.1)  as 


D0y/Dt  =  3rr  w  +  Z 


(5.7) 


where  for  brevity  in  what  follows  we  let  3=0g» 

The  ensemble  average  equations  in  the  second-order  closure  system 
invol ving  0^  through  Equation  (5.7)  are  those  governing  0^,  0^0^,  and 
q'0y.  These  moments,  derived  from  Equation  (5.7)  and  the  corresponding 
conservation  laws  for  the  other  variables  are 


D0v/Dt  =  -  3(u!0;)/axi  +  3r(r'w'  +  r  w)  +  Z 


(5.8) 


+  (gi/Tr)ev'^-2eijk«jU,;0/  +  0. 75(q/A)u ‘0^ 


+  0.3  9(qA3u!0^/9Xj)/aXj  +  3rj  ujujr 
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(5.9) 


+  3rj  Uj  ujr'  +  E'U] 


D(e;2)/Dt  =  -2uj0;(30y/9Xj-35r)-O.45(5/A)0;2 


+  0.3  9(qA30;2/axj)/9Xj 


+  3rj  UjOyr'  +  3rj  uj  0^r'  +  E'0^ 


(5.10) 


D(q'0;)/Dt  =  -  ujq'  (9?y/9xj-3rj7)  -  uj0;  9q/9xj 
-  0.45(q/A)  q'0;  +  3rj  ujq'r' 

+  0.3  9(qA9q'9^/9xj)/9Xj  +  3rj  Uj  ^'r‘  +  i~q '  (5.11) 

In  the  above  equations  =  r  gj/  g  ,  and  q  and  A  are  characteristic  velocity 
and  length  scales  for  the  turbulence. 

We  observe  that  the  effective  buoyancy  driving  the  production  terms  is 
the  generalized  buoyancy  9^  defined  as 

Qw  =  0v  ■  PrR  (5.12) 

where  R  is  defined  as  the  cloud  depth 

z 

R(z)  =  j"  r(z)dz  (5.13) 

Zo 

and  z  is  the  coordinate  direction  aligned  with  the  gravitational  body  force. 
From  Eqs.  (5.12)  and  (5.13)  we  then  have 
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(5.14) 


ae^^/3xj  =  ae^/axj-arj-r 


D9^/Dt  =  Dey/Dt-3rrw  (5.15) 

Equation  (5.7)  expressed  in  terms  of  the  generalized  buoyancy  thus  becomes 
conservative' (except  for  the  radiation  and  precipitation  source  term  E): 


De^^/Dt  =  E 


(5.16) 


The  generalized  buoyancy  0^^  is  an  appropriate  conserved  variable  for 
partly  cloudy  situations.  Although  0^  is  conserved,  it  is  important  to  note 
that  0y  (not  remains  as  the  determinant  buoyancy  for  the  momentum 
equation. 


It  will  be  observed  that  new  second  and  third  order  correlations 
involving  r  are  introduced  into  the  system.  After  the  mean  cloudiness  r,  the 
most  important  second  order  correlation  is  the  cloud  flux  w'r' .  The  cloud 
flux  appears  directly  as  a  heating  source  term  in  the  mean  virtual  potential 
temperature  Equation  and  carries  the  essence  of  cumulus  heating  in  partly 
cloudy  situations.  A  third  order  correlation,  particularly  important  in 
conditionally  unstable  flows,  is  the  cloud  vertical  turbulence  energy  1/2 
w'w'r'  which  represents  the  vertical  turbulence  energy  which  is  correlated 
with  the  cloud  fluctuations.  Similarly  important  are  the 
cloud  heat  and  moisture  fluxes  w'0y'r',  w'q'r'  which  represent  the  fluxes  of 
heat  and  moisture  which  are  embedded  in  the  fluctuating  clouds.  (Note  that 
all  these  correlations  vanish  in  uniformly  saturated  stratus  for  which  r  =  1, 
r'  5  0.) 


The  ensemble  mean  liquid  water  qj^  and  variance  (which  do  not  appear 
directly  in  Eqs.  (5.8)-(5.11)  but  control  the  radiative  transport  and 
precipitation  microphysics)  are,  assuming  liquid-vapor  equilibrium. 


(5.17a) 
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+  2  (q^  q'r'2  +  r 


(5.17b) 


It  can  be  seen  that  the  liquid  water  variance  (a  second-order  correlation) 
involves  third  and  fourth  order  correlations  of  the  basic  variable  q^=q-q5. 

Since  r  itself  can  be  directly  expressed  in  terms  of  q^,  it  is  not  a  new 
dependent  variable  of  the  flow;  rather  it  is  a  function  of  the  basic  set  (u^, 

'l)*  Thus,  the  above  moments  can  be  expressed  in  terms  of  moments 
involving  only  the  basic  set.  If  the  distribution  function  of  the  variables 
over  the  ensemble  were  known,  this  would  then  be  a  simple  matter  of  formal 
calculation.  We  shall  show  how  the  partly  cloudy  correlations  can  be  reduced 
to  expressions  involving  only  the  fundamental  set  and  certain  coefficients 
which  relate  cross-correlations  and  variances  of  the  partial  cloudiness  r  with 
the  fundamental  variables.  To  proceed  further  to  explicit  forms  for  these 
coefficients,  further  information  about  the  distribution  functions  must  be 
specified.  We  have  carried  through  the  complete  explicit  evaluations  of  the 
coefficients  and  all  correlations  for  the  case  of  Gaussian  distributions;  we 
shall  quote  the  results  of  the  evaluations  below. 

We  now  consider  the  reduction  of  the  new  correlations  involving  r 
appearing  in  Eqs.  (5.8)-(5.11)  to  expressions  in  terms  of  the  fundamental  set 
of  variables  (u^-,  0y,  q).  For  equilibrium  systems,  the  quantity  q;^  =  q-q^  is 
central  to^  these  partly  cloudy  correlations.  From  Oliver,  Lewellen,  and 
Williamson  ,  the  fluctuation  q^  may  be  expressed  in  terms  of  the  fundamental 
set  as 


q{  =  aq'-be; 


(5.18) 


where 
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(5.19a) 


3  =  1-  yq^ms 

b  =  YgBjqs/Tr  (5.19b) 

and 

Yq  =  -  0.61  +  (y"^+0.61)r  (5.20a) 

Ye  =  1  +  G~^-'^)T  (5.20b) 

The  correlations  and  (j)'q^  can  now  be  expressed  in  terms  of  the  fundamental 
variables  as 


'^qx  "  =  a2  q'2  -  2ab  q'e^  +  b2  0^2 

(5.21) 

i>'q{  =  a  q'4,'  -  b  0;<j,' 

(5.22) 

Let  us  define  the  coefficients  by  the  statement 

=  C(^)  (op/cq^)  (5.23) 

where  cp  =  (r'r')^/2  is  the  cloudiness  variance.  Equation  (5.23)  is,  in 
effect,  the  defining  Equation  for  the  coefficient  For  variables  <j),  q^ 

which  are  Gaussian,  the  coefficients  are  unity. 

We  may  then  express  the  partly  cloudy  second-order  correlations  in  terms 
of  the  fundamental  variables  (using  Eqs.  (5.22)-(5.23))  as 

ujr'  =  c('^i)  (®p/<Jq^)  (a  u^  -  b  (5.24a) 

e;r'  =  c(®)  (cJr/%)  (a  ^)  (5.24b) 
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=  C(q)  (or/oq  )  (a  q'2  -  b  q'e;) 

A 


(5.24c) 


q'r' 


The  correlation  Oq^  is  given  by  Equation  (5.21)  in  terms  of  the 
fundamental  set.  The  representation  of  Op  and  the  coefficients  require 
further  information  about  the  distribution  function.  We  summarize  the  results 
for  Gaussian  variables  here.  The  cloudiness  state  Q  is  defined  as 

Q  =  qx/%  (5.25) 

and  represents  the  ratio  of  the  mean  saturation  difference  to  the  rms 
saturation  difference.  For  Q-x»  the  state  must  approach  one  of  full  cloudiness 
(r-»-l)  and  vanishing  fluctuations  (op^O).  For  Q-»— »  the  state  must  be  fully 

clear  air  (r-»-0)  with  similarly  vanishing  fluctuations  (ap-»-0).  When  Q=0  the 

mean  flow  is  just  at  the  saturation  point  and  the  fluctuations  should  be  at  a 
maximum  (op'^^rmax)*  Gaussian  distributions  of  any  variable  ^  and  the 

difference  mixing  ratio  q^,  it  can  be  shown  that 

7  =  [1  +  Erf  (Q//2)]  (5.26) 


/2) 


exp  (-Q2/2) 


(5.27) 


e  1  (5.28) 

Hence  the  cloud  state  Q  is  determined  by  the  fundamental  mean  and  second-order 
variables;  and  this  single  parameter  in  turn  determines  Op,  r,  and  The 

cloudiness  variables  are  thus  closed  in  terms  of  the  fundamental  variables. 

From  Equation  (5.24a)  the  cloud  flux  is  represented  in  terms  of  the  total 
moisture  flux  and  heat  flux  as 

w'r'  =  (op/oq^)  (a  w'q'  -  b  w'e^)  (5.29) 
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The  cloud  flux  represents  the  net  flux  resulting  from  partly  cloudy  updrafts 
and  down  drafts.  When  w'q'  and  w'9y  have  the  same  sense  the  cloud  flu?<  is 
diminished  by  the  presence  of  heat  flux  due  to  warming  which  is  positively 
correlated  with  moistening.  When  w'Qy  opposite  sense  to  w'q'  the  cloud 
flux  is  enhanced  by  the  positive  correlation  of  cooling  with  moistening. 

Third-order  qorrelalfions  appearing  in  Eqs.  (5.9)-(5.11)  present  more 
tasks  in  modeling.  Following  Donaldson  (1973),  we  may  possibly  model  these 
terms  as  gradient-driven  diffusion  of  second-order  correlations: 

u^w'r'  -  Vp  (qA)  au^r'/az  (5.30a) 

OyW^r'  *=  -  Vp  (qA)  a  Q'r'/az  (5.30b) 

q'w'r'  T  -  Vp  (qA)  a  q'r'/3z  (5.30c) 

The  diffusion  coefficient  Vp  may  be  taken  as  equal  to  v^.  the  general  diffusion 
coefficient  in  the  second-orfler  closure  system.  However,  we  also  recognize 
that  this  simple  approach  may  not  be  adequate  for  intermittent  turbulent 
flows. 


Turbulence  Closure  Theory  on  an  Intermittent  Turbulence  Field 

The  modeling  procedure  and  modeling  coefficients  developed  for 
second-order  closure  thqory  rest  on  the  assumption  that  the  smallest 
resolvable  regions  of  the  flow  field  under  consideration  are  fully  turbulent 
(or  fully  laminar).  In  large  scale  meteorological  models,  a  conditionally 
unstable  atmospheric  flow  is  often  in  a  quite  different  state.  It  may  be 
composed  mostly  of  stable  regions  of  vanishingly  small  turbulence  which 
co-exist  with  a  small  fraction  of  highly  turbulent  zones.  We  believe  that 
second-'Order  modeling  can  quite  accurately  describe  the  kinds  of  turbulent 
features  present  in  these  turbulent  zones  such  as  plumes  and  jets,  penetrating 
downdrafts,  and  gust  fronts  if  the  detailed  structure  of  these  features  is 
resolved.  On  the  other  hand^  if  we  choose  to  resolve  only  the  large  scale 
flow  which  consists  at  any  instant  of  a  small  fraction  of  turbulent  zones 
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embedded  in  fluid  which  is  otherwise  weakly  turbulent,  we  must  describe  such  a 
field  by  introducing  an  intermittency  function  There  are  a  variety  of 

ways  of  defining  the  intermittency;  however,  we  shall  use  a  simple  partition 
definition.  Let  q  denote  a  volume  of  flow  which  is  small  compared  to  the 
large  scale  gradients  of  interest  yet  contains  a  number  of  disturbed  zones 
within  it.  Let  denote  the  fraction  of  Q  which  contains  the  disturbed 
regions,  while  denotes  the  undisturbed  portion.  The  intermittency 

is  then  defined  as 


(5.31) 

In  this  preliminary  investigation  we  shall  temporarily  set  aside  the  mean 
velocity  field  effects  and  illustrate  what  amounts  to  the  free  convection 
limit.  The  mean  velocity  and  corresponding  shearing  production  effects  can 
later  be  incorporated  with  the  same  methodology. 

The  governing  equations  of  the  large  scale  variables  are  the  equations 
for  the  virtual  potential  temperature  and  water  mixing  ratio  with  the  partly 
cloudy  correlations  presented  in  Part  5. 


D9v/Dt  =  -  9w'e;/3z  +  Br(rw+r'w' )+E 

(5.32) 

D’q/Dt  =  -  9w'q79z-C 

(5.33) 

In  the  above  C  is  the  mean  precipitation  rate  and  E  contains  the  effects  of 
radiative  transport  as  well  as  precipitation. 

Because  of  intermittency,  the  second-order  closure  system  described  in 
Part  2  must  be  modified  to  describe  the  large  scale  correlations  w'S^,  w'r', 
w'q' .  We  assume,  however,  that  in  the  disturbed  regions  of  fraction  o)2 , 
turbulent  correlations  may  be  defined  which  are  governed  by  the  usual 
second-order  closure  equations.  Let  (  denote  an  average  for  the 

disturbed  regions.'  The  usual  second-order  equations,  assumed  descriptive  of 
these  disturbed  regions,-  are  then  of  the  typical  form 
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D(w'e;)d/Dt  =  -  (w'w')d  (8ew/3z)d  +  9/Tp(e;2)^  .  0.75(qH/A)(w'9 ■ )h 


+  0.3  3qHA8(w'9  '  )^/3z+3r(w'w' r '+w  w '  r '  ' w '  )(j  (5.34) 

Since  we  assume  the  turbulence  is  negligible  in  the  fraction  (l-w^)  of 
the  flow  field,  the  large  scale  higher-order  correlations  (  )  may  be 

related  to  the  disturbed  region  correlations  (  as 

(  )  =  a)2(  (5.35) 

and  for  the  rms  quantities, 

q  =  coq^  (5.36a) 

Op  =  a)(ap)jj  (5.36b) 

X  "  (5.36c) 

The  large  scale  mean  cloudiness  is  given  by 

r  =  a)2(7)(j+(l-a)2)(7)g  (5.37) 


where  (r)g  is  the  environmental  region  cloudiness  of  fraction  (l-a)2). 

Let  us  now  summarize  the  developments  required  to  apply  the  second-order 
closure  system  to  partly  cloudy  and  possibly  intermittent  flows.  The  first 
development  concerns  the  dissipation  and  conventional  triple  correlation 
(diffusion)  terms  which  are  indicated  by  a  single  underline  in 
Equation  (5.34).  We  expect  that  the  characteristic  velocity  q  in  these 
modeled  terms  should  be  based  on  the  average  turbulent  velocity  in  the 
turbulent  region  rather  than  on  the  average  over  the  larger  scale  region. 
Consistent  with  Equation  (5.36a),  it  therefore  appears  that  one  of  the 
principal  modifications  required  in  the  system  of  equations  to  make  them 
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dpplicsble  to  dn  intormittont  flow  is  to  divido  tho  siriQlG  undorliriGd  tGrms  by 
0)  so  that  these  modeled  terms  are  now  based  on  the  appropriate  characteristic 
velocity.  Thus,  the  large  scale  second-order  correlation  equations  in 
intermittent  flow  may  be  expected  to  be  modeled  as 

Dq^/Dt  =  g/Tr^  -  a)"l  (q/4A)q^+a)-l0.33qA8qV82)/3Z  (5.38) 


D(e(,2)/Dt  =  -  w'e^  (90y/3z-3rr)+  0.3a)"l  3qA80y2/3z 

-  a)"iO.45q/A)(0^2)+3p(^  0;,r '+w'0^,r'  )  +  (z  *0^.)  (5.39) 


Dw'q'/Dt  =  -  w'w'3q/3z+(g/Tp)  q'0^+O.3co"l  3qA3w'q73z) 

-  0.75a)"l(q/A)w'q'  (5.40) 


Dw'0;/Dt  =  -  w'w'  39^/3z  +  (g/Tp)0;2  _  0.75  l-l(q/A)w'0; 

+  0.3  0)"!  3(qA3w'0'/az)3z  +  3r(w'w'r'  +  w  w ' r ' )  +  i  'w'  (5.41) 


D(q'e;)/Dt  =  -  (w'q')((8ev/92)-3rr)  -  (w'9;)(3q/3z) 

+  O.3a)"l3qA3q'0^/3z)/3z 

-  0.45aj“i  (q/A)q'ey  +  3rw  q  V  +w '  q  '  r '  +  Z  'q  '  (5.42) 


D(q'2)/Dt  =  -(w'w')d(3q/9z)d-0-45a)"l(5/A)q'2 
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+  0.3(jj"l3(qA8q'2/3z)/3z 


(5.43) 


New  modeling  may  be  required  for  the  doubly  underlined  terms  which  are 
triple  correlations  involving  cloud  fluctuations.  In  strongly  intermittent 
flow,  the  gradient-driven  modeling  previously  given  may  well  be  inadequate  for 
these  terms.  It  is  possible  that  in  such  flows  a  model  of  the  form 


w'(t)'r'  =  apOj'^Op  w'<|) 


(5.44) 


(where  ap  is  an  order  unity  modeling  constant)  may  be  appropriate.  Such  a 
form  has  the  correct  limiting  behavior  in  that  it  vanishes  in  clear  air  or 
stratus  situations  when  ap->0  and  in  strongly  intermittent  flow  it  has  the 
property 


w'(|)'  r'  ~  w'(() 


since  ((jp)jj~l  and  thus  uj'^Op  ~  1  in  strongly  intermittent  flow,  as  indicated 

by  Equation  (5.36b).  Another  possibility  is  that  some  more  global  integral 
type  model  may  be  more  appropriate  for  these  terms. 

The  third  development  required  is  a  theory  of  closure  for  the 

intermittency  itself.  In  our  discussion  up  to  this  point  we  have  utilized 

the  intermittency  simply  to  distinguish  the  disturbed  regions  from  the 

2 

undisturbed  regions.  We  now  present  the  conditions  which  determine  w  .  For 
absolutely  unstable  flows  in  which 


3e^/3z  <  0 


(5.45) 


the  flow  is  non-intermittent  and  we  set  u=l  for  such  flows. 


For  flows  which  are  conditionally  stable,  we  have 


0  <  39^^/3z  <  (l-r)3r 


(5.46a) 


r  <  1 


(5.46b) 
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For  conditionally  unstable  flow  we  offer  the  following  tentative  closure  for 
the  intermittency.  Convection  is  induced  by  the  condensation  of  moisture. 
Thus,  the  convectively  active  zones  are  necessarily  the  cloudy  zones.  As  the 
mean  cloudiness  declines,  a  point  is  reached  in  which  the  partly  cloudy 
updrafts  and  downdrafts  cannot  fill  the  whole  space  and  a  bifurcation  into 
disturbed  and  undisturbed  regions  occurs. 

In  the  conditionally  unstable  regime,  the  flow  will  be  unstable  to 
vanishingly  small  disturbances  for  all  >  0.  At  =  0,  the  stability  is 
neutral.  We  postulate  that  this  neutral  stability  point  is  also  the  point  of 
bifurcation  into  disturbed  and  undisturbed  regions.  For  q;^  <  0  and  30^^/3z 
satisfying  Equation  (5.46a)  the  flow  is  intermittent  with  (q;^)g  <  0  in  the 
undisturbed  regions.  We  postulate  that  in  this  intermittent  regime  (^x)d 
remains  at  the  level  of  the  flow  at  the  bifurcation  point;  thus 

(qx)d  “lx  °  (5.47) 

That  is,  the  disturbed  or  cloudy  region  should  be  characterized  by  a  mean 
moisture  level  which  is  approximately  equal  to  the  local  saturation  level. 
From  Equation  5.26  it  then  follows  that 

I'd  =  "2  ’  ""  ~Jy  (5.48) 

consistent  with  identifying  the  cloudy  regions  with  the  disturbed  regions  we 
must  take  the  undisturbed  zones  as  cloud  free.  Then  from  Eqs.  (5.36),  (5.37) 
and  (5.48) 


r  =  03^2,  0.=-^  (5.49) 

/2 

Closure  may  now  be  completed  by  recognizing  that  some  fraction  of  the 
humidity  fluctuations  in  the  disturbed  region  must  be  of  sufficient  magnitude 
to  exceed  the  amount  by  which  the  environmental  humidity  departs  from  the 
saturated  area.  Thus 
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(5.50) 


where  Cq  now  becomes  a  modeling  coefficient  determined  by  what  fraction  of  the 
fluctuations  are  assumed  to  exceed  -(qx)  •  If  this  is  set  at  the  5%  level 
then  our  assumption  of  a  Gaussian  distribution  in  the  disturbed  region  would 
give  Cq  =  0.6.  As  long  as  Cq  is  assumed,  Eqs.  (5.36c)  and  (5.50)  then 
determine  the  intermittency 


(5.51) 


where  we  have  used  (q^)  =  qx/(l-“  )• 


e 


The  system  of  equations  is  now  complete  with  the  mean  variables 
determined  by  Eqs.  (5.32)  and  (5.33).  The  second-order  correlations  of  the 
mean  variables  are  determined  by  Eqs.  5.38  to  5.43;  the  cloud  flux  term  w'r' 
by  Equation  5.29;  r  and  by  Equation  5.49;  by  Equation  5.21  and 

finally  (u  by  Equation  5.51. 

For  application  to  the  problem  of  sub-grid  cumulus  description  (cumulus 
parameterization),  we  would  expect  to  develop  simplified  approximate  solutions 
to  this  set.  Some  of  the  simplifications  which  may  be  possible  include  the 
approximations  of  quasi-equil ibrium,  super-equilibrium,  or  self-similarity. 
The  most  appropriate  approximations  for  cumulus  parameterization  remain  to  be 
determined.  It  is  important  to  note  that  the  only  inputs  to  the  turbulence 
equations  are  the  large  scale  mean  fields  and  the  input  fluxes  at  the  cloud 
base  (or  the  full  boundary  layer  may  be  solved  concurrently).  From  these 
inputs  the  cloud  flux  w'r'  should  emerge  naturally  yielding  the  heating  rate 
di stribution. 
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Data  Comparisons 


To  examine  the  consistency  of  this  modeling  procedure,  we  make  some  order 
of  magnitude  analyses  with  an  actual  data  set.  The  data  set  chosen  is  that 

3  9 

for  the  Gate  B  Scale  Array  (Thompson  et  al .  ).  Selected  observations  from 

this  data  set  are  tabulated  in  Table  1. 


TABLE  1.  SELECTED  DATA  FROM  THE  GATE  B  SCALE  ARRAY 
(Thompson  et  al ) 

(w'q’)o 

(from  precipitation)  1.4  x  lO'**  m/s 

h 

region  of  conditional  stability  10*+  m 

(950  mb  to  185  mb) 

aq/az 

(characteristic  moisture  -1.8  x  10"6  m"i 

gradient  in  mid  portion 
of  layer  (~h/2)) 

90'v/3z 

(characteristic  buoyancy  +  0.0035  °C/m 

gradient  in  mid  portion 
of  layer  (~h/2)) 

T*”  ^ 

Low  Level  Convergence 

and  Evaporation  Rate  7.5  x  10"^  s"^ 

gr 

(in  mid  portion  of  0.0040  °C/m 

layer  ~h/2) 

qx 

-0.002 

Average  observed  heating  rate 

by  cumulus  convection  5°C/day 

Although  the  cloudiness  r  is  not  measured,  it  is  observed  to  be  small. 
From  Table  1  it  can  therefore  be  seen  that  the  system  is  conditionally  stable, 
since  gr  >  ae^/az  >  0.  Since  the  latent  heat  release  is  important  in 
conditionally  stable  systems,  the  buoyancy  flux  is  determined  principally  by 
moisture  flux.  Let  (w'q')o  be  the  characteristic  moisture  flux.  The 
characteristic  buoyancy  is  then  (w'ey)Q  ~  (L/Cp)(w'q' )q.  From  the  turbulence 
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kinetic  energy  Equation  (5,38),  we  have  in  super-equilibrium 


a)-lq3/A  ~  (g/Tp)  w'e^; 


(5.52) 


hence  we  may  construct  a  characteristic  velocity  scale  w+  as 


w3  =  (g/Tp)(L/Cp)(w'q')o  a)h 


(5.53) 


where  h  is  the  thicknqss  of  the  conditionally  stable  region. 

The  characteristic  scale  of  the  moisture  fluctuation,  q+,  is  defined  as 


q+  =  (w'q')o/w+ 


(5.54) 


Then  from  Equation  5.51 


(1)  =  - 


<1+ 


O-H  0-6qx 


(5.55) 


73 


0) 


(w'q')o 


2/3 


“•H  (f  ^ 

\  r  ‘'p 


73 


(5.56) 


For  the  data  values  given  in  Table  1,  we  find 


j-1/3  =  4.9  m/s 


ji/3  =  0.29xl0”‘*  m/s 


0)  =  0.06 


w^  =  1.9  m/s 
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q+  =  0.7x10"'+ 


r  =  1.9x10"3 


=  0.024 

The  value  of  mean  cloudiness  for  this  data  set  seems  a  bit  low;  however,  the 
cloudiness  variance  is  of  the  order  of  2%. 

If  the  cumulus  flux  is  approximated  from  the  leading  term  of 
Equation  (5.29)  as 


w'  r ' 


a(ar/0q^) 


(w'q')o 

2 


With  a  ~  1  and  the  above  values,  we  obtain  w'r'  =  0.02  and  the  cumulus  heating 
rate  is  then 


erw'r'  =  8°C/day 

which  is  as  close  to  the  observed  heating  rate  as  we  should  expect  for  the 
current  rough  approximations. 

As  an  independent  check  of  the  intermittency,  we  may  calculate  the 
moisture-flux  in  mid  layer  from  super  equilibrium  using  the  observed  moisture 
gradient  and  compare  it  with  the  precipitation  derived  value  (w'q ' )q.  In 
super-equil ibrium 


I  I  ■  *  A 

w'q'  ~  -ojqA  8q/9z 

with  A=103m,  q  ~  w^  ~  1.9,  and  9q/9z  ~2xl0”®  m"^  we  find 

w'q'  ~  1.8x10"'+  m/s 
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which  corresponds  with  the  observed  1.4x10”**  m/s  surface  value. 

Although  we  have  not  yet  used  this  set  in  a  predictive  model,  this 
exercise  does  indicate,  at  least,  an  order  of  magnitude  consistency  for  this 
new  set  of  equations. 

We  believe  that  this  is  a  promising  approach  to  the  problem  of  cumulus 
parameterization.  This  is  particularly  true  for  mesoscale  models,  for  which 
the  Arakawa  and  Schubert  scheme,  most  favored  for  global  models,  would  have 
difficulty  in  regions  where  the  intermittency  is  not  much  less  than  one. 
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6.  CONCLUDING  REMARKS 


We  believe  that  the  results  of  the  Cal  span  fog  simulations  are  as 
accurate  as  should  be  expected  for  a  one-dimensional  simulation  of  what  may 
generally  be  expected  to  be  a  three-dimensional  phenomenon.  Further, 
substantial  improvements  in  modeling  the  atmospheric  boundary  layer  in  the 
marine  environment  appear  to  require  addressing  directly  the  problem  of  three 
dimensionality.  This  is  particularly  true  of  the  coastal  regions  which  are 
vital  to  the  Navy's  interests.  We  have  the  choice  of  either  developing  a 
complete  unsteady,  three-dimensional  version  of  our  second-order  closure 
model,  or  using  the  existing  one-  and  two-dimensional  versions  of  our  model  to 
develop  appropriate  turbulent  flux  parameterization  schemes  which  may  be  used 
to  adequately  specify  the  sub-grid  processes  in  three-dimensional,  regional 
meteorological  models  developed  by  others.  As  discussed  in  Chapter  IV,  we 
believe  the  latter  approach  is  preferable  at  the  present  time,  because  it 

should  permit  our  model  developments  to  have  an  earlier  impact  on  operational 
models. 

The  grid  resolution  in  any  mesoscale  meteorological  model  will  never  be 
adequate  to  completely  resolve  the  turbulent  transport  processes  of  importance 
in  the  troposphere.  Thus,  accuracy  of  the  physical  mechanisms  which  are 
controlled  by  turbulence  will  depend  on  how  faithfully  the  sub-grid 

parameterization  can  simulate  these  processes.  Our  second-order  closure  model 

of  turbulent  transport  in  the  atmosphere  provides  two  ways  of  developing  this 
parameterization.  First,  there  are  the  dynamic  equations  for  the  second-order 
flux  quantities  of  interest,  which  may  be  approximated  in  some  fashion. 

Second,  our  existing  one-  and  two-dimensional  models  of  the  atmospheric 
boundary  layer  may  be  run  with  relatively  fine  resolution  to  test  the  accuracy 
of  any  tentative  parameterization  schemes. 

In  order  to  investigate  how  well  a  candidate  parameterization  scheme  may 
be  integrated  into  a  mesoscale  model,  we  expect  to  use  the  three-dimensional 
Navi er-Stokes  model  which  Dr.  Sykes  has  currently  running  on  A.R.A.P.  s 
computer.  This  model  currently  exists  in  two  forms:  as  a  model  of  the  full 
Navier-Stokes  equations  in  cartesian  coordinates  as  used  by  Mason  and  Sykes 
and  as  a  hydrostatic  model  in  a  terrain-following  coordinate  system  as 
proposed  by  Clark**^-  This  model  should  serve  as  a  test  vehicle  to  determine 
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the  relative  merits  of  increasingly  complex  parameterization  schemes. 

Some  of  the  most  important  turbulent  transport  processes  in  the  lower 

troposphere  involve  clouds.  Therefore,  in  order  for  any  sub-grid  flux 

parameterization  scheme  to  be  very  successful  on  the  mesoscale  level,  it  must 

incorporate  some  of  the  essential  elements  of  cloud  dynamics.  As  discussed  in 

Chapter  V,  v/e  have  reviewed  the  cumulus  parameterization  schemes  proposed  by 
„  ,  36  38 

Arakawa  and  Schubert  and  Kuo  .  It  appears  that  the  Arakawa  and  Schubert 
model  which  has  found  some  success  in  global  models  would  require  extensive 
modifications  to  be  used  at  the  mesoscale  level.  On  the  other  hand,  we 
believe  our  second-order  closure  approach  can  provide  a  more  realistic 
representation  than  the  Kuo  model  which  is  much  simpler  than  that  of  Arakawa 
and  Schubert. 
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APPENDIX  A 

MODELING  THE  ROLE  OF  TURBULENCE  IN  CLOUD  MICROPHYSICS 

INTRODUCTION 

In  this  work  we  present  a  model  of  cloud  and  warm  precipitation  which  is 
naturally  coupled  to  the  turbulence  of  the  atmosphere  in  which  clouds  reside. 
Models  of  precipitation  previously  developed  have  been  constructed  for 
specific  cloud  and  storm  systems  (Kessler,  1969;  Wilhelmson  and  Klemp,  1981; 
Orville  and  Chen,  1982)  and  take  no  mechanistic  account  of  the  effects  of 
turbulence.  In  the  present  work  we  explicitly  describe  certain  of  the  actions 
of  turbulence  upon  the  growth  processes  of  cloud  droplets.  These  include  the 
effects  of  turbulence-induced  shear  and  acceleration  fields  on  the  droplet 
coagulation. 

In  Part  2  we  review  the  stages  of  condensation,  cloud  droplet  spectrum 
evolution,  and  precipitation;  and  we  point  out  the  role  of  turbulence  in  the 
various  microphysical  processes.  In  Part  3  we  present  a  microphysical  closure 
model  which  reduces  the  integro-di fferential  kinetic  equation  describing  the 
droplet  population  to  a  tractable  form  which  explicitly  incorporates  the 
effects  of  turbulence  described  in  Part  2.  In  Part  4  we  develop  the  turbulent 
ensemble  average  forms  of  the  cloud  and  precipitation  equations  and  the 
second-order  correlation  variables  according  to  the  methodology  of  higher 
order  closure  theory.  A  particular  development  of  the  ensemble  average  is 
required  for  description  of  the  turbulent  motion  of  the  cloud  and 
precipitation  drops  in  the  turbulent  cascade  and  the  corresponding  ensemble 
average  of  droplet  functions.  In  Part  5  we  illustrate  some  of  the  properties 
of  the  model  for  homogeneous  clouds. 
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CLOUD  FORMATION,  EVOLUTION,  AND  PRECIPITATION 


To  set  the  stage  for  the  microphysical  closure  we  shall  develop,  we  first 
review  the  various  processes  which  take  place  from  the  onset  of  a  water  mixing 
ratio  in  excess  of  the  saturation  value  to  the  final  stage  (if  it  occurs  in 
the  time  scale  of  a  particular  macro-event)  in  which  drops  precipitate  to  the 
surface.  We  are  particularly  interested  in  the  role  of  turbulence  in  these 
various  stages  of  development.  The  stages  of  drop  evolution  may  be  defined  in 
the  scheme  of  Table  Al. 

Once  the  water  mixing  ratio  exceeds  the  local  saturation  value,  nuclei 
must  be  activated  before  water  may  condense  in  realistic  time  scales. 
Following  nuclei  activation,  drops  grow  by  the  direct  condensation  of  vapor 
and  are  col  1 isional ly  coupled  only  through  Brownian  motion.  In  most 
atmospheric  situations  the  liquid  water  formed  by  the  overall  amount  of  excess 
mixing  ratio  over  the  saturation  value  and  the  number  of  nuclei  available  and 
activated  results  in  a  cloud  with  drop  number  densities  ranging  from 
10^-10^  m“3  and  radii  ranging  from  1  to  10  ym. 

A  particular  feature  of  the  growth  of  cloud  droplets  and  the  evolution  of 
the  spectrum  is  the  creation  of  a  small  number  of  drops  much  larger  than  the 
average.  This  tail  effect  in  the  distribution  function  is  intensified  by  the 
large  collision  cross-section  of  large  drops,  the  mechanism  by  which  droplets 
grow  significantly  beyond  the  range  of  radii  <  lOym  up  to  precipitation  sizes 
in  excess  of  100  ym  is  still  an  outstanding  unresolved  problem.  One  class  of 
mechanism  is  collisional.  Since  Brownian  collisional  rates  are  much  too  small 
for  significant  collisions  over  characteristic  macro  time  scales  for  this  size 
range  of  droplet,  gravitational  sedimentation  has  long  been  identified  as  a 
collisional  coalescence  mechanism  of  atmospheric  clouds,  (Berry,  1967; 
Warshaw,  1967;  Twomey,  1966).  There  seems  little  doubt  that  the  collisional 
coalescence  of  droplets  of  different  size  by  gravitational  sedimentation  is  an 
important  droplet  growth  mechanism  at  some  stage  of  cloud  evolution. 
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TABLE  Al.  STAGES  OF  CLOUD  DROPLET  SPECTRUM  FORMATION  AND  EVOLUTION 


(1)  The  Nuclei  Activation  Stage 

drop  (nucleus)  size 
response  time  scale 

(2)  The  Condensation  Growth  Stage 

drop  size 

response  time  scale 

(3)  The  Collisional  Growth  Stage 

drop  size 

time  scale 


10“^  Mm  -  1  urn 
1-10  sec 


1  pm  -  20  pm 

10  sec 


10  ym  -  10"^  pm 

2  4 

10  -  10  sec 


(4)  The  Sedimentation  Stage 
drop  size 
time  scale 


50  pm-10“3  ym 
0.1  -  1  hr 


On  the  other  hand,  in  the  early  stages  of  growth  (radii  from  5-30  pm)  this 
mechanism  possesses  certain  limitations.  Two  of  these  limitations  are  the 
inherent  requirement  of  differential  size  for  a  non-zero  collision  rate  and 
the  sharply  diminished  collision  efficiencies  which  result  for  the  low 
Reynolds  numbers  of  droplet  sedimentation  in  the  range  of  sizes  from  5  to 
30  pm  (Mason,  1971).  These  two  limitations  when  viewed  in  the  light  of  a 
further  result  from  classical  condensation  theory  -  namely  the  narrowing  of 
the  droplet  spectrum  during  condensation  growth  (Sedunov,  1974;  Levin  and 
Sedunov,  1966)  -  suggest  that  growth  mechanisms  other  than  gravitational 
sedimentation  may  play  an  important  role  in  the  growth  of  cloud  droplets  into 
precipitation  size  drops. 

It  would  appear  that  atmospheric  turbulence  can  play  an  important  and 
direct  role  in  the  evolution  of  the  cloud  droplet  spectrum  through  both 
collisional  and  non-col  1 i sional  processes.  The  collisional  processes  include 
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the  coagulation  induced  by  turbulent  shearing  and  acceleration  fields  (Saffman 
and  Turner,  1955).  The  non-coll isional  processes  involve  modification  of  the 
droplet  spectrum  by  either  fluctuations  or  turbulent  mixing  of  the  humidity 
and  temperature  during  condensation  and  evaporation.  This  latter  effect  has 
been  studied  in  a  statistical  manner  (Levin  and  Sedunov,  1966;  Sedunov,  1974; 
Jeou  Jong,  1966)  and  in  simplified  homogeneous  cloud  models  (Bartlett  and 
Jonas,  1971;  Mason  and  Jonas,  1974;  Baker,  et  al.,  1980;  Jonas  and 
Mason,  1982).  These  studies,  which  are  not  without  contradiction  of  one 
another,  indicate  that  turbulent  fluctuations  of  supersaturation  in  absence  of 
mixing  do  not  necessarily  broaden  the  cloud  droplet  spectrum.  On  the  other 
hand,  turbulent  entrainment  of  environmental  air  into  cloud  coupled  with  large 
scale  growth  and  decay  fluctuations  on  the  scale  of  the  lifetime  of  individual 
clouds  do  appear  to  appreciably  broaden  the  droplet  spectrum.  Such 
non-col  1 isional  "turbulent"  broadening  of  the  spectrum  would  clearly  augment 
the  coll isional  processes  which  ultimately  complete  the  evolution  of  cloud 
into  precipitation  size  drops. 

The  final  stages  of  drop  evolution  occur  when  drops  have  grown  large 
enough  to  develop  a  significant  sedimentation  velocity.  These  precipitable 
drops  then  leave  the  cloud  and  progress  to  the  surface  where  they  leave  the 
atmosphere.  Turbulence  may  often  play  a  significant  role  at  this  stage.  We 
shall  show  that  under  most  atmospheric  conditions,  drops  as  large  as  1000  urn 
are  still  strongly  correlated  with  the  largest  turbulent  eddies.  Turbulent 
contributions  to  the  sedimentation  flux  of  precipitation  must  therefore  be 
included  in  the  construction  of  the  total  flux  of  water  drops  which  reach  the 
surface.  In  general,  the  turbulent  flux  can  either  oppose  or  augment  the 
sedimentation  flux  and  lead  to  a  diminishment  or  enhancement  of  the  net  flux. 

On  the  basis  of  the  time  scales  presented  in  Table  Al,  we  may  regard  the 
small  size  spectrum  dynamics  (nuclei  activation  and  condensation  -  evaporation 
rates  in  the  elementary  spectrum)  as  occurring  instantaneously  on  the  time 
scale  of  turbulence.  Correspondingly,  the  collisional  growth  and 
precipitation  evaporation  time  constants  are  often  much  greater  than  the 
turbulence  time  scale.  We  shall  use  this  difference  in  time  scales  in  the 
construction  of  an  approximate  description  of  the  microphysical  processes  and 
in  simplifying  the  turbulent  ensemble  equations. 
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DESCRIPTION  OF  CLOUD  AND  PRECIPITATION 


Cloud  and  Precipitation  Moment  Equations 

Instead  of  considering  the  detail  of  the  entire  droplet  spectrum,  we 
divide  the  spectrum  into  two  groups:  cl oud,  consisting  of  small  droplets 
which  do  not  precipitate,  and  precipitation,  consisting  of  large  droplets 
which  possess  a  sedimentation  flux. 

The  choice  of  the  spectrum  dividing  volume  VpQ  which  separates  cloud 
drops  from  precipitation  must  be  made  on  the  basis  of  two  principal 

constraints.  The  first  is  that  VpQ  lie  in  the  tail  region  of  the  elementary 

spectrum  (i.e.,  the  spectrum  before  collisional  processes  become  active)  thus 
we  require  VpQ>a^VQ,  where  Vq  is  the  average  elementary  volume  and  is  the 
dispersion  of  the  cloud  elementary  spectrum.  Second,  since  the  small  droplets 
have  negligible  sedimentation  velocities  (over  macro  time  scales)  we  do  not 
specify  the  precipitation  flux  as  an  average  over  the  entire  liquid  water 

distribution.  Rather,  it  is  more  appropriate  to  define  the  precipitation 

group  as  those  droplets  with  sedimentation  velocities  greater  than  a  certain 
minimum  value.  This  minimum  value  is  determined  by  the  overall 
macro-dynamics.  This  sedimentation  velocity  is  selected  so  that  a  drop  will 
fall  over  a  characteristic  macro-length  in  some  characteristic  macro-time. 
For  example,  for  cloud  with  a  length  scale  of  the  order  of  1  km  or  less  and 
precipitation  on  a  time  scale  of  the  order  of  1  hour  or  less,  the  minimum 
precipitation  velocity,  VpQ  should  be  of  the  order  of  10-100  m/hr  which 
corresponds  to  the  sedimentation  velocity  of  a  droplet  of  approximately 
20-30  ym  radius. 

Since  the  elementary  drop  size  lies  in  the  range  of  4-10  ym  for  most 
clouds,  it  appears  that  the  rule  of  thumb  Vpp/Vg  ~  10-20  may  be  adequate  for 
most  situations  yielding  a  spectrum  dividing  volume  outside  the  bulk  of  the 
elementary  cloud  distribution,  yet  not  being  so  great  that  a  significant 
portion  of  the  smallest  precipitation  drops  are  inadequately  represented  nor 
that  the  bulk  of  the  cloud  spectrum  becomes  included  in  the  precipitation 
portion  of  the  spectrum. 
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The  number,  mass,  and  size  relationships  within  the  cloud  and 
precipitation  are  then  as  follows.  Let  Vq  be  the  smallest  droplet  volume  of 
interest.  Let  m=(v/vQ)  be  a  size  specification  parameter  based  upon  droplet 
volume  and  hence  proportional  to  liquid  mass.  The  spectrum  dividing  size 
parameter  is  mpQ=(VpQ/vQ).  We  denote  the  number  density  of  cloud  drops  of 
size  m  as  n^-Cm)  and  of  precipitation  drops  as  np(m).  The  cloud  and 
precipitation  distribution  functions  f(.(m),  fp(m)  are  then  defined  as 

fj-(m)  =  h(.(m)/n(.  l<m<mpQ  (A.l) 

fp(m)  =  np(m)/n^  mpQ<m  (A. 2) 

where  n^  and  np  are  the  total  number  densities  of  cloud  and  precipitation, 
respectively  (Figure  Al).  The  distribution  functions  thus  satisfy  the 
normalization  conditions 

f(.  (m)  =  1  (A. 3a) 


fp  (m)  =  1  (A. 3b) 

The  mixing  ratio  of  cloud  and  precipitation  qp  are  given  by 

qc  =  (Po/P“)''o^c"’c  (A. 4a) 

qp  =  (Po/P“)''oPp"’p  (A. 4b) 

where  the  liquid  water  density  is  Pq,  the  density  of  the  air-cloud  mixture  is 
p„  and  m^.,  mp  are  the  average  sizes  of  cloud  and  precipitation: 
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Figure  A.l  Schematic  of  liquid  drop  distribution  and  cloud  and 
precipitation  groups.  Source  droplets  exist  at  size  m=l.  The 
cloud  group  spans  the  range  1  <  m  <  m  Bi -modality,  when  it 
occurs,  occurs  with  the  second  mode  in  the  precipitation  group. 
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The  cloud  and  precipitation  dispersion  are  defined  as 


a 


2 

c 


(rn^-nic)  (m) 


(A. 5b) 


(A. 6a) 


2  2  2  2 

2^  (m  -trip)  fp  (m)  (A. 6b) 

"’="ipo 

Let  q  be  the  total  water  mixing  ratio  and  the  mixing  ratio  of  vapor. 
The  mixing  ratios  satisfy  the  conservation  condition 

q  =  qy  +  q^-  +  qp  .  (A. 7) 

Utilizing  the  fact  that  the  nuclei  activation  and  condensation  growth  rates  of 
droplets  of  sizes  less  than  50  ym  are  so  rapid  and  that  cloud  will  generally 
fall  well  within  this  range,  we  assume  that  cloud  is  in  liquid-vapor 
equilibrium  over  the  macro-time  scale.  Since  precipitation  may  involve  large 
drops  with  slow  evaporation  time  scales  in  cloud-free  air,  we  shall  not 
require  the  precipitation  to  be  in  equilibrium.  The  presence  of  cloud  will 
thus  turn  on  the  level  of  the  saturation  mixing  ratio  q^  relative  to  the  total 
mixing  ratio  q.  We  introduce  the  cloudiness  r(qx)  which  is  a  conditional 
variable  indicating  the  presence  or  absence  of  cloud.  It  is  defined  as  the 
Heaviside  function  of  the  difference  mixing  ratio  qx=q-qs  as 
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r  =  H(qx)  . 


(A.8) 


The  equilibrium  liquid  water  mixing  ratio  q^  is  then  given  by 


qt  =  Hx  •  (A. 9) 

The  cloud  and  vapor  mixing  ratios  q^.,  q^  are  given  in  terms  of  q,  qp,  q^  as 

qc  =  (qx-qp)  .  (a.io) 

qv  =  q-^qx-(i-^)qp  •  (a.h) 

We  shall  take  as  the  dynamical  equations  governing  cloud  and  precipitation  the 
first  two  moments  of  the  droplet  kinetic  equation  for  each  group:  the  number 
density  and  mixing  ratio  moments.  (We  use  the  dynamical  equation  for  q  rather 
than  q^  since  either  serves  as  an  appropriate  dependent  variable).  The 
results  are 


Dq/Dt  +  9(qpVpi)/8Xi  =  0  ,  (A. 12) 

Dn^/Dt  -  “N(;(;-Np^-N^p+N^y  ,  (A. 13) 

Dqp/Dt  +  9(qpVpi)/axi  =  M^p+Mpc-Mp,  ,  (A.14) 

Dnp/Dt  +  3(npVp^.)/9Xi  =  N^p-Npp  .  (A.15) 


In  the  above  N(-(.+Nj.p  represents  the  loss  of  cloud  droplets  due  to 
self-collisions  among  cloud  droplets;  N^-p  and  M^-p  are  the  production  of 
precipitation  number  and  mixing  ratio  by  self-collisions  within  cloud 
(auto-conversion  of  cloud).  Np^.  and  Mp^.  are  the  rate  of  loss  of  cloud  number 
and  the  rate  of  production  of  precipitation  mixing  ratio  by  collisions  of 
precipitation  drops  with  cloud  drops  (cloud  collection).  Npp  is  the  rate  of 
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loss  of  precipitation  number  density  by  self  collision  of  precipitation.  For 
heavy  precipitation  with  drops  in  excess  of  2  mm,  collisional  de-coalescence 
and  droplet  instability  break-up  terms  must  be  included  in  Equations  A. 13 
through  A. 15.  and  Mp^  are  droplet-vapor  interaction  terms. 

represents  the  net  production  of  cloud  number  due  to  condensation  and 

evaporation.  Mp^  represents  the  evaporative  rate  of  precipitation  in 

cloud-free  air. 

The  variables  q,  n^,  qp,  Op  and  Eqs.  A. 12  through  A. 15  represent  the 
basic  dynamical  variables  and  conservation  laws  for  cloud  and  precipitation  in 
our  simplified  description.  The  cloud  and  vapor  mixing  ratios  are  determined 
from  Equations  A. 10  and  A. 11  and  the  sizes  m^.  and  mp  from  Equation  A. 4.  The 
precipitation  speed  then  follows  from  the  sedimentation  velocity  function 
V^(m)  as 


(A. 16) 


and  the  precipitation  flux  as  qpVp^- . 

The  closure  of  Equations  A. 12  through  A. 15  requires  specification  of  the 
various  collisional  rates  N,  M  and  the  liquid-vapor  rates  Mpy.  The  exact 
representation  of  the  collisional  rates  N^-p,  and  M^-p  are 


N 


cc 


K(m,k)fc(m)fc(k) 


(A. 17a) 


N 


cp 


K(m,k)fj.{m)f(.(k) 


(A. 17b) 


%  =  (nipo/mj  (q^/nj  N^p  (A.17c) 

where  K(m,k)  is  the  collision  kernel  between  a  droplet  of  size  m  and  a  droplet 
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of  size  k.  We  may  express  these  results  in  terms  of  a  non-dimensional 
collisional  effectiveness 


N 


cc 


n^K 


C'CC  ""CC 


(A. 18a) 


^cp  "  '^c^cp  ^cp  (A. 18b) 

where  K^-j-  and  K(-p  are  appropriately  defined  average  collision  kernels.  The 
collisional  effectiveness  is  a  measure  of  the  effectiveness  of  al  1  sizes 
within  cloud  in  evolving  the  spectrum  to  precipitation  size  drops; 


S 


cc 


K(m,k)  f(-(m)  f(-(k) 


(A. 19a) 


m„^-2  tn_Q-l 


^cp  "  *^cp 


■  2 


K(m,k)  f^(m)  f^(k) 


m=l  k=mpQ-m 


We  may  similarly  express  the  rates  Nr,r>  and  Nr,-,  as 

r'-  PP 


^pc  '^p'^c  *^pc  ^pc  ’ 


(A. 19b) 


(A. 20a) 


^pc  ”  '^c’^p  ^pc  ^pc 


(A. 20b) 


Npp  =  Hp  Kpp  Spp 


(A. 20c) 


where  Sp^.  and  Spp  are  appropriate  collisional  effectiveness  for  the  collection 
of  cloud  by  precipitation  and  for  precipitation  self-collisions. 
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The  collision  kernel  K(m,k)  appearing  in  the  collision  rate  expressions 
is  composed  of  shearing  and  differential  sedimentation  parts.  Since  we  do  not 
describe  the  detailed  microphysics  of  the  small  size  spectrum  formation  but 
only  the  collisional  growth  stage  for  droplets  greater  than  1  urn  radii,  we 
neglect  Brownian  collisions.  The  collision  kernel  is  thus  expressed  as 

K(m,k)  =  K(^)(m,k)  +  K(^^)(m,k)  ,  (A. 21) 

where  K(^)(m,k)  is  the  kernel  for  droplets  in  a  shear  field  and  is  the 
kernel  for  differential  sedimentation  in  a  force  or  acceleration  field.  These 
kernels  have  the  form 

K(s)(m,k)  =  Vo/Tr(m^/3+k^/3)3  e  e(5)(m,k)  .  {A.22) 


K(ds)(m.k)  =  7r(3Vo/4TT)^/3  (m'/3+k'/3)"(i-p„/pQ)  Vi(m)-Vi(k)  e(ds)(m.k)  , 

(A. 23) 

where  G  is  the  magnitude  of  the  fluid  shear,  is  the  droplet  drift  velocity, 
and  e(^)(m,k)  and  e(^^)(m,k)  are  the  collisional  efficiencies. 

The  efficiencies  for  differential  sedimentation  are  given  by 

Mason,  1971.  The  collisional  efficiencies  in  shear  flow  are  less  well 
established.  For  particles  of  nearly  identical  size  e^^^  -  1.  Swift  and 
Friedlander,  1964,  suggest  a  value  e(^)  -  0.4. 

The  sedimentation  velocity  function  V^(m)  may  be  expressed  in  terms  of 
Tp(m)  as 


Vi(m)  =  (l-P„/Po)  Tp(m)  +  91^  ,  (A. 24) 

where  Du^/Dt  is  the  fluid  acceleration,  g^-  the  acceleration  of  gravity,  and 
Tp(m)  is  the  droplet  relaxation  time. 
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For  the  purposes  of  the  present  theory,  the  droplet  relaxation  time  in 
air  at  the  earth's  surface  may  be  adequately  approximated  as 


T^(m)  =  aR^ 


(A. 25) 


where  is  the  radius  of  a  drop  in  size  class  m  and  a,  a  are  constants  given 
in  Table  A2  which  is  constructed  from  the  data  of  Gunn  and  Kinzer  (1949): 


TABLE  A2 

Droplet  Radius  Rj^ 

_ urn _  _ a _  _ ^ 

1  <  Rm  <  40  1.6x10^  s/m^  2 

40  <  R^  <  600  6.4x10^  s/m  1 

600  <  R^  <  1500  20  s/m^/2 


Selected  values  of  are  given  in  Table  A3. 


TABLE  A3 


R 

(ym) 

(I 

10 

1 

50 

32 

100 

64 

500 

320 

1000 

632 

Microphysical  Closure 

The  elementary  volume  Vq  as  well  as  the  collisional  and  liquid-vapor 
rates  in  Equations -A. 12  through  A. 15  involve  the  full  distribution  functions 
f(,(m),  fp(m)  consistent  with  the  fact  that  a  reduced  description  of  the  full 
spectrum  to  two  groups  consisting  of  cloud  and  precipitation  cannot  be  a 
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closed  representation  without  additional  closure  approximation.  Three  overall 
properties  of  the  full  spectrum  are  critical  for  the  macroscopic  cloud  and 
precipitation  dynamics.  These  are  the  elementary  volume  Vq,  the  dispersion  of 
the  cloud  spectrum  a^,  and  the  structure  of  the  cloud  distribution 

1  in  the  vicinity  m  mpg.  The  elementary  volume  Vq  establishes  the 
average  size  of  the  cloud  drop.  The  cloud  dispersion  g  as  well  as  the 
average  size,  Vq  are  the  determining  microphysical  properties  for  the 
collisional  evolution  of  the  bulk  of  the  cloud  spectrum  and  hence  are  critical 
to  the  rate  N^.^.  The  structure  of  the  tail  of  the  cloud  distribution  function 
determines  the  population  of  tail  droplets  and  hence  is  critical  to  the  rates 
Ncp.  M(-p  or  correspondingly  the  effectiveness  S^-p. 

Rigorous  closure  of  the  elementary  volume  requires  a  dynamical 
description  of  both  the  condensation  nuclei  and  the  supersaturation.  In  lieu 
of  such  a  detailed  description,  we  fix  the  elementary  volume  by  specifying  a 
parameter:  the  elementary  number  density  n^,  which  may  be  thought  of  as  the 

effective  population  of  cloud  droplets  in  the  absence  of  collisions.  The 
elementary  volume  is  then  determined  by  ng  as 


(A. 26) 


which  from  Equation  A. 4a  then  determines  the  cloud  size  m^.  as 


"ic  =  V^c  •  (A. 27) 

The  elementary  number  density  may  be  thought  of  as  primarily  an 
environmental  and  cloud  type  parameter.  Typical  values  are  given  in  Table  A4. 
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TABLE  A4 


Cl oud/Envi ronment  Elementary  Number  Density  (cm~3) 

Continental  Cumulus  200-800 

Maritime  Cumulus  100-200 

Maritime  Stratus  50-100 

The  presence  of  collisions  which  we  shall  discuss  in  what  follows  acts  as 
a  broadening  mechanism  for  the  spectrum  and  leads  to  an  increasing  dispersion 
a^.  The  manner  in  which  non-col  1 isional  processes  increase  (or  decrease) 
is  not  presently  clearly  understood.  We  shall  therefore  restrict  attention  to 
the  case  where  the  cloud  spectrum  in  absence  of  collisions  consists  of 
droplets  concentrated  at  a  single  size  m(.=l  and  a  negligible  dispersion  0. 
Once  collisions  come  into  play,  we  have  m^-  >  1  and  >  0. 

Since  the  condensation  and  evaporation  rates  of  the  elementary  spectrum 
are  rapid  compared  to  macroscopic  time  scales,  we  assume  that  the  elementary 
drop  number  appears  instantly  when  the  saturation  line  is  crossed. 
Correspondingly,  the  rate  N^-y  may  be  represented  as 

where  <S(q^)  is  the  Dirac  delta  function.  In  cloud  free  air  the  rate  Mpy  is 
given  by 

^pv  ~  ^p/^pv  (A. 29) 

with  the  rate  Xpy  given  by 

—  [l+CRe/2Sc'/3]  (q^-q^)  (A.30) 

Sc 

where  n  is  the  kinematic  viscosity  of  air,  and  Rg,  are  Reynolds  number  and 
Schmidt  number  based  upon  droplet  diameter  2Rp  and  sedimentation  velocity  Vp. 
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0:20.28  is  a  correlation  constant.  Since  qv=qs  when  r=l  it  can  be  seen  that 
Mpy=0  except  in  clear  air. 

We  now  consider  the  col  1 i sional ly  generated  tail  of  the  distribution 
function  which  will  allow  us  explicit  closure  for  and  S(-p.  The 

collisional  processes  will  also  determine  It  can  be  shown  that  even  for  a 

constant  collision  kernel,  the  effectiveness  S^-p  is  a  sensitive  function  of 
the  tail  of  the  distribution.  This  can  be  demonstrated  for  the  special  case 
of  a  Smoluckowski  cloud  (droplet  population  initially  n(.(0)  of  a  single  size 
m=l  (O(-=0)  at  time  t=0  under  the  action  of  a  constant  collision  kernel  K 
(Smoluckowski,  1917). 

For  such  a  cloud,  the  distribution  function  is  given  by 

f^  (m,T)  =  fc  (1,T)T"’"1  (A. 31) 

where 

fc  (1,T)  =  (1-T)/(1-t"’p°'S  (A. 32) 

and  the  function  T  for  the  Smol uchowski ' s  problem  is  given  by 

T  =  (l/2)n^(0)Kt/[l+(l/2)n^(0)Kt]  (A. 33) 

Upon  substitution  of  Equation  A. 26  into  the  definition  of  A. 19,  we  find 

Sec  =  '"c^c(l»T)  (A. 34a) 

^cp  =  ( V"’c)^c("’po-l»''')  (A. 34b) 

The  effectiveness  S^.^.  which  describes  cloud-cloud  collisions  which  do  not 
create  precipitation  and  S^-p  which  describes  cloud-cloud  collisions  which 
create  a  precipitation  drop  are  shown  in  Figure  A2  as  a  function  of  the 
evolution  parameter  T.  It  can  be  seen  that  at  T=0,  we  have  S^p=0  since  no 
droplets  which  are  concentrated  at  f^(l,0)=l  can  reach  the  spectrum  dividing 
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Figure  A. 2  Properties  of  the  cloud  distribution  function  as  a  function  of 
the  evolution  parameter  T. 


size  in  a  single  collision.  Correspondingly  at  T=0,  we  have  since  all 
droplets  which  exist  at  f(,(l,0)=l  will  remain  as  cloud  droplets  after  their 
first  collision.  At  T=l,  we  have  a  flat  distribution  with  f(-(m,l )=l/(mpQ-l ) 
and  with  S(.(.=S^p=l/2-l/(mpQ-l ) .  For  T  >  1,  the  distribution  shifts  with  most 
cloud  droplets  near  the  spectrum  dividing  size  where  f(;(mpQ-l)  ~  1.  In  this 
state  we  have  S^.^.  0  and  S^-p  -  1. 

Description  of  the  evolution  of  the  distribution  function  is  thus 
essential  for  the  proper  description  of  cloud  evolution  and  the  precipitation 
production  rate.  As  our  collision  rate  closure,  we  shall  assume  that  the 
cloud  distribution  function  to  be  used  in  the  determination  of  is  a  power 
law  form  consistent  with  the  Smoluchowski  form  of  Equations  A. 31  and  A. 32.  We 
regard  the  function  T  as  an  evolution  function  which  describes  the  evolution 
of  the  cloud  droplet  spectrum  from  one  which  is  peaked  about  the  elementary 
spectrum  size,  v=Vq,  m=l,  T=0  to  one  which  is  flat  with  all  sizes  equally 
represented:  f^(m)  l/(mpQ-l).  We  fix  the  value  of  T  through  Equation  A. 27. 
For  the  form  Equation  A. 31,  m^(T)  from  Equation  A. 5a  is  given  by 

mc(T)  =  [l-(mpo-l)fc(l)  T>-1]  (A.35) 


Thus,  the  evolution  function  T  is  determined  implicitly  by  Equations  A. 27  and 
A.35.  The  dispersion  for  the  form  Equation  A. 31  is  given  by 


2,  . 


“r.(T)  ■ 


po 

ni^(l-T) 


(A. 36) 


Under  these  closure  hypotheses  the  collisional  rates  become 


^cc  '^c'^cc^cc 


(A. 37a) 


%  =  (A.37b) 

^cp  “  '^c^c'^cc^cp(^)  (A. 37c) 
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(A.37d) 


^pc  ~  ‘^p^c'^pc 
^pc  ~  *^p^c^pc  (A.37e) 

^pp  “  "2  '^p^pp  (A.37f) 

Because  of  the  strongly  diminished  collision  efficiency  e^^^  for  particles  of 
significantly  disparate  size,  we  neglect  the  contributions  of  shear  collisions 
to  Kpj..  Because  of  the  large  differential  sedimentation  rate  of  large  drops, 
we  similarly  neglect  shear  contributions  to  Kpp. 

The  effectiveness  with  variable  collision  kernel  is  assumed  to  be 
approximately  the  same  as  the  effectiveness  with  a  constant  collision  kernel. 
This  closure  assumption  places  some  burden  upon  an  appropriate  choice  for  the 
average  kernels  Kp^,  Kpp.  We  select  as  "average"  kernels  the  forms 


K^.(.  =  K^^)(mJ,m")  +  K(^^)(mJ,m")  , 

(A. 38a) 

Kpc  =  K(^s)(mp,mJ  , 

(A. 38b) 

Kpp  =  K(^s)(m+,mp  , 

(A. 38c) 

where  the  sizes  are  specified  as 

m*  =  m^(l±a^./2)  , 

(A. 39a) 

""p  =  • 

(A. 39b) 
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The  precipitation  dispersion  Op  is  not  a  critical  parameter  of  the  model 
since  it  principally  controls  the  rate  Npp.  The  rate  Npp  is  not  a  critical 
rate  since  the  precipitation  flux  divergence  9(Vp^qp)/9x^  will  dominate  over 
Npp  for  precipitation  drops  in  excess  of  50  urn  in  clouds  of  the  order  of 
several  km  or  less  in  depth.  We  assume  Op  is  given  in  terms  of  the 
precipitation  flux  by  a  parameterization  (Marshall  and  Palmer,  1948);  it  may 
be  as  satisfactory  to  take  it  as  a  fixed  cloud  type  parameter.  We  assume  that 
the  precipitation  distribution  function  shapes  are  not  important  for  the  rates 
Npc,  Npp  and  it  is  permissible  to  set  Sp(.=l,  Spp=l/2. 

Since  the  power  law  form  Equation  A. 31  is  always  monotone  decreasing  with 
size  it  does  not  allow  bi -modality  within  cloud.  It  should  be  noted,  however, 
that  the  collisional  process  which  induces  bi-modality  into  the  spectrum 
(other  than  specially  chosen  initial  conditions)  is  the  rapid  increase  of  the 
collision  kernel  (both  in  collision  cross-section  and  efficiency)  with  an 
increase  in  size  of  the  larger  collision  partner.  There  is  thus  a  consistency 
between  the  cut-off  of  the  monotone  decrease  of  the  cloud  spectrum  at  m=mpQ 
and  the  large  collision  kernel  associated  with  drops  for  m>mpQ.  Hence  the 
total  droplet  spectrum  allows  for  bi-modality  with  the  second  mode  (when  it 
occurs)  always  occurring  in  the  precipitation  part  of  the  spectrum 
(Figure  Al). 

Our  microphysical  theory  is  thus  closed.  The  moisture  variables  in  a 
cloudy,  precipitating  atmosphere  are  described  by  the  moisture,  cloud  and 
precipitation  variables  q,  q„,  n_,  n-,.  These  variables  are  governed  by  the 

pup 

conservation  laws  Equations  A. 12  through  A. 15.  The  rates  appearing  in  these 
conservation  laws  are  given  by  Equations  A. 37.  The  evolution  function  T  is 
central  to  these  collisional  and  liquid-vapor  rates  and  is  determined  by  m^ 
from  Equations  A. 31  and  A. 35.  The  model  coefficient  of  the  closure  is  the 
elementary  number  density  ng. 
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TURBULENT  ENSEMBLE  AVERAGE  EQUATIONS 
Mean  Equations 

The  turbulent  ensemble  averages  of  the  cloud  and  precipitation  equations. 
Equations  A. 12  through  A. 15,  are 

Dq/Dt  +  aCiJT^  +  Vp^.qp  +  V'^.q')/9Xi  =  0  ,  (A. 40) 


Dn^/Dt  +  3(Uin,)/3x,.  =  -Ncc^N^p-Np^+Ncv 


{A.41) 


Dqp/Dt  +  3  (u'q-+Vp^.qp+V^q')/3Xi  = 


pv 


(A. 42) 


Dnp/Dt  +  3(u>‘+Vpinp+V'in')/3xi  =  N^p-Npp  .  (A. 43) 

Since  Vp^-=V^(mp)  and  mp  is  a  function  of  qp  and  Op  through  Equation  A. 25, 
the  fluctuation  Vp^-  may  be  represented  as 

''pi  =  ^1%  -  »  (A.44) 


where 


3V  • 

-l£l 

3V  ■ 
b.  =  —21 

3mp 

»qp  ■ 

1  3mp 

3np 

(A. 45) 


The  correlations  involving  Vp  then  reduce  to  correlations  involving  the 
fundamental  set  of  variables  q,  qp,  n^-,  np  and  the  turbulent  precipitation 
fluxes  of  mixing  ratio  and  number  may  be  expressed  as 


V 


-  bi 


qp^p 


(A. 46a) 


V 


qp^p 


bi 


(A. 46b) 
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For  Vp^  ~  it  follows  from  Equations  A. 46a  and  A. 46b  that 


^pi*lp  ~  (^p  /^p  "  ^p'^p/'lp^p)  ^pi*lp  »  (A. 47a) 

_  _ j. _ 

^pi'^p  ~  (^p^p/*lp*^p  ”  ’^p  /*^p)  ^pi  *^p  •  (A. 47b) 

We  shall  term  the  flux  Vp^qp+Vp|qp  the  precipitation  flux  and  the  flux 
u^qp  the  turbulent  precipitable  water  flux.  We  see  that  net  flux  consists  of 
both  turbulent  and  precipitation  contributions.  The  term  u-jqp  represents  the 
turbulent  transport  of  precipitation  drops  by  turbulent  updrafts  and  down 

draftsj _ U  can  significantly  diminish  or  enhance  the  pure  precipitation  flux 

'^pi*lp''''^pi‘lp* 

The  new  non-col  1 i sional  correlations  introduced  which  directly  involve 
cloud  and  precipitation  are  the  turbulent  precipitable  water  flux  u]qp,  the 
turbulent  cloud  number  flux  u-jn^,  and  the  turbulent  precipitation  number  flux 
uiOp.  The  collisional  correlations  required  are  Np^.  M^p,  and 

Npp.  The  condensation-evaporation  correlations  required  are  and  Mpy. 

The  ensemble  averages  of  the  derived  moisture  variables  follow  from 
Equations  A. 10  and  A. 11  as: 

^c  =  ^  -  ^'^p  >  (A.48a) 

qy  =  r  qs  +  (l-r)(q-qp)  -  r'q;^  -  r'qjl,  ,  (A. 48b) 

^X  =  q  -  ^  •  (A. 48c) 

We  see  that  there  are  second-order  correlation  contributions  to  the  mean  cloud 
and  vapor  mixing  ratios. 
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First  Order  Rate  Approximation 


On  the  basis  of  the  slow  time  scales  of  the  collisional  rate  terms  in 
Eqs.  A. 13  through  A. 15,  we  propose  a  first  order  rate  approximation  in  which 
we  retain  the  mean  collision  rates  N(-p,  Np^-,  M^-p,  Mp^.,  Npp  and  the 

liquid-vapor  rate  Mpy  (which  must  also  appear  in  the  virtual  potential 
temperature  equation),  but  we  neglect  their  fluctuation  contributions  to  the 
second-order  correlation  equations  for  u]qp,  u'n^L,  u^Op,  etc. 

In  this  approximation,  the  instantaneous  rates  given  by  Equations  A. 37 
become 


^CC  “  ^CC  » 

(A. 49a) 

^cp  ~  ^CC  ^cp  ’ 

(A. 49b) 

^cp  ’’  ^'^c^c^'^c^c^  ^cc  ^cp  ’ 

(A. 49c) 

^pc  "  ^pc  * 

(A.49d) 

^pc  "*  (^^p^C^'^p^C^  ^pc  * 

(A.49e) 

^pp  ~  2  ^ 


PP 


(A.49f) 


In  the  above  decomposition  we  have  assumed  the  droplet  loading  is  never  large 
enough  to  affect  the  eddy  dynamics  and  thus  the  collision  kernels  are 
de-correlated  from  the  droplet  concentrations.  We  have  also  interpreted  S^.^. 
as  (T),  etc.  The  ensemble  average  of  the  collision  kernels  given  by 

Equations  A. 22  through  A. 25  require  the  ensemble  average  fluid  shear  G  and  the 
ensemble  average  net  acceleration 


g-j+Du-j/Dt  .  For  simplicity,  but  without 
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essential  restriction,  we  shall  neglect  the  contribution  of  mean  fluid  shear 
3u^/3Xj  and  mean  acceleration  Du^/Dt.  Following  Saffman  and  Turner  (1955),  we 
approximate  the  ensemble  average  acceleration  as 


gi+Du)/Dt  =  /g2+a2  (A. 50) 

where  a  is  the  rms  acceleration 


Before  we  can  evaluate  the  quantities  a  and  G  as  well  as  all  correlations 
involving  droplet  variables  n^.,  qp,  Op  we  must  first  consider  the  structure  of 
atmospheric  turbulence  and  its  influence  on  the  motion  of  liquid  drops. 
Atmospheric  turbulence  consists  of  the  motion  of  eddy  structures  ranging  from 
the  largest  energy  containing  scale  A  to  the  microscale  where  molecular 
dissipation  comes  into  play.  Correspondingly,  the  large  scale  eddies  have  a 
characteristic  time  scale  Ty^=A/q  where  q  is  the  rms  turbulence  velocity. 
The  turbulence  dissipation  rate  e=l/4q^/A  is  preserved  through  this  turbulent 
cascade.  The  smallest  eddies  of  scale  XQ=(n^/e)  possess  a  characteristic 
microtime  =(ri/E)  a  characteristic  shear  -  z  1'^  and  a 

,  0  3/00 

characteristic  acceleration  a^  =  -  e  The  magnitudes  of  these 
quantities  for  the  range  of  turbulence  dissipation  rates  encountered  in  the 
atmosphere  are  presented  in  Table  A5. 
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TABLE  A5.  TURBULENCE  LENGTH,  TIME,  ACCELERATION,  AND  SHEAR  SCALES 
IN  THE  ATMOSPHERE  FOR  DISSIPATION  SCALE  EDDIES 


e 

^0 

^Xo/9 

^X 

'^O 

(m^/sec^ ) 

(ym) 

(ms) 

- 

(s'M 

0.001 

1510 

130 

0.009 

7.7 

0.01 

846 

42 

0.05 

23.8 

0.10 

476 

13 

0.28 

76.9 

1.0 

268 

4.2 

1.58 

238 

10.0 

151 

1.3 

8.9 

769 

Since  the  dissipation  rate  is  preserved  through  the  cascade, 
fluctuation  time  of  an  eddy  of  scale  x  is  (x/Xq) 

characteristic  shear  in  an  eddy  of  scale  X  is 
acceleration  in  an  eddy  of  scale  X  is  a^=ap^  (X/X  )" 


the 

The 

the 


For  those  conditions  where  Tp/Tj^>>l  the  droplet  motion  de-correlates 
completely  from  the  turbulence  of  the  air  and  there  is  no  influence  of 
turbulence  on  the  collision  kernels  appearing  in  Equations  A. 22  and  A. 23.  The 
ensemble  mean  shear  G-  and  acceleration  a  experienced  by  the  drops  in  the 
absence  of  mean  shear  and  acceleration  tends  to  zero.  At  the  opposite  extreme 
for  «1  the  droplets  are  completely  correlated  with  the  full  turbulent 

cascade.  The  corresponding  ensemble  mean  shear  G  which  such  droplets  will 
experience  is  given  by 


G  =  /2/15  Gx  .  (A. 52) 

0 

(Taylor,  1935).  The  rms  acceleration  which  the  droplet  will  experience  is 


a  =  /I. 3  a^ 


(A. 53) 


(Batchelor,  1950). 


131 


Let  us  now  examine  the  droplet  relaxation  time  relative  to  the  eddy 
fluctuation  time  From  Table  A3  it  can  be  seen  that  over  a  considerable 

range  of  importance  the  droplet  relaxation  time  is  much  greater  than  the 
microscale  time  and  the  motion  of  such  droplets  will  not  be  correlated 
with  the  motion  of  dissipation  scale  eddies.  On  the  other  hand,  the  macrotime 
is  rarely  less  than  10  seconds  under  atmospheric  conditions  (except  deep  in 
the  surface  layer)  so  that  the  droplet  motion  is  correlated  (even  for  1000  ^m 
drops)  with  a  large  part  of  the  turbulent  cascade. 

When  we  postulate  that  the  smallest  eddy  scale  X  whose  motion 

the  droplet  will  follow  is  determined  by  the  condition  and  thus  this 

scale  is  given  by 

^  =  ^0  (A. 54) 

correspondingly  the  maximum  shear  Gx  and  maximum  acceleration  which  such 
drops  can  experience  is 


Gx  =  GXq  (VAXq)'^  .  (A.55) 

^X  =  ^Xq  •  (A.56) 

We  now  formulate  the  general  ensemble  average  shear  and  accelerations  as 


G(m)  =  c(^)  xTl  Min 
'^0 


(A. 57) 


a(m)  = 


a^  Mi  n 


(A. 58) 


where  we  indicate  the  explicit  dependence  upon  droplet  size  and  c^^^  and 
are  model  constants.  For  droplets  which  are  embedded  in  an  isotropic  cascade 
we  may  take  c(^)=  /2/15  ,  c(‘^^)=/l,3. 
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It  can  be  seen  that  for  Tp(m)^  ,  the  ensemble  mean  shear  G(m)  becomes 

independent  of  e;  the  ensemble  mean  acceleration,  however,  continues  to 

increase  with  e  as  a(m)  ^  a^  whereas  a(m)  a^  for 

.  .  '"o  ^0  ^0 

Tp(m)  ^  . 

0 

The  ensemble  average  kernels  are  thus  given  by  Equations  A. 38  and 

Equations  A. 22  through  A. 24  with  G  given  by  G(m^)  and  Du^/Dt  by  a(m^)  where  m^ 

is  the  size  of  the  larger  collision  partner  and  G(m),  a(m)  are  given  by 
Equations  A. 57  and  A. 58. 

The  relative  importance  of  turbulent  shearing  collisions  compared  to 

gravitational  sedimentation  collisions  in  evolving  the  elementary  spectrum 
depends  upon  the  magnitudes  of  the  turbulence  dissipation  level  e  and  the 
droplet  spectrum  dispersion  In  Figure  A3  we  exhibit  the  ratio  of  the 

turbulent  shear  kernel  to  differential  gravitational  sedimentation  kernel 

In  Figure  A4  we  show  the  ratio  of  the  total  collisional  kernel 

with  all  turbulence  effects  included  to  the  turbulence-free 
differential  sedimentation  rate  K^^s),  ^an  be  seen  in  Figure  A3  that  for 

cloud  droplets  dispersed  over  the  range  of  6  to  20  ym  radius,  turbulence 
induced  shearing  collisions  play  an  important  role  in  cloud  evolution  for 

dissipation  rates  in  excess  of  0.01  m^/s^.  For  wider  dispersions, 
differential  sedimentation  will  be  more  dominant;  for  narrower  distributions 
turbulence  shearing  dominates  the  collisional  process.  For  strong  turbulence 
as  might  exist  in  cumulus  corresponding  to  e  ''  1-10  m^/s^,  turbulence  induced 
shearing  collisions  may  be  the  principal  collisional  broadening  and  spectral 
evaluation  mechanism.  In  Figure  A4  it  can  be  seen  that  when  both  turbulent 
shearing  and  acceleration  effects  are  included  in  the  full  collision  kernel, 
turbulence  effects  dominate  the  collision  kernel  for  e  i  0.05  m^/sec^. 

Turbulence  Equations 

In  the  first  order  rate  approximation,  we  neglect  the  contributions  of 

the  collisional  and  liquid  vapor  rates  to  the  higher-order  turbulence 

correlations.  The  turbulence  equations  for  the  cloud  and  precipitation 
variables  u]n(I,  u^np,  u]qp  and  their  other  appropriate  second-order 
correlations  are  those  of  second-order  closure  theory  for  passive  scalars. 
(Lewellen,  1978;  Oliver,  Lewellen  and  Williamson,  1978).  Since  drops  up  to 
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Figure  A. 3  Isopleths  of  the  ratio  of  the  turbulent  shear  collision  kernel 

1^0  the  differential  sedimentation  kernel  (in  which 
the  turbulence  induced  acceleration  a  is  suppressed)  as  a 
function  of  turbulence  dissipation  rate  e  and  the  difference  in 
collisional  partner  radii  AR/R  .  The  elementary  drop  radius  R„ 
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1000  ym  are  still  correlated  with  the  large  eddies,  the  turbulence  equation 
for  the  cloud  and  precipitation  will  not  require  any  de-correlation  of  the 
droplet  motion  in  the  bulk  of  the  flow.  Appropriate  turbulence  equations  with 
droplet  de-correlation  will  only  be  required  deep  in  the  surface  layer  where 
the  time  of  turbulence  approaches  the  relaxation  time  of  the  largest  drops. 
Such  surface  layer  forms  have  been  developed  by  Lewellen  (1977). 

Homogeneous  Cloud  Illustration 

The  full  properties  of  the  present  theoretical  model  for  cloud  and 
precipitation  can  only  be  revealed  when  it  is  integrated  with  a  fully-coupled 
inhomogeneous  turbulence  model  such  as  that  of  Oliver,  Lewellen,  and 
Williamson  (1978).  We  can,  however,  reveal  certain  of  its  features  including 
the  general  behavior  of  the  evolution  time  and  precipitation  intensity  in 
simple  homogeneous  cloud  illustrations.  In  this  first  homogeneous  cloud 
illustration,  we  neglect  all  transport  terms  in  Equations  A. 40  through  A. 43 
except  the  precipitation  flux  divergence  which  we  approximate  as 

with  a  similar  representation  for  ^(Vp-jOp+Vp^  Op)  /3x^  where  is  a 

characteristic  cloud  depth.  We  further  assume  that  at  time  t=0  the  total 
liquid  water  existing  is  cloud  water  so  that  qp(0)=0.  We  also  assume  that  the 
liquid  water  existing  initially  is  not  replenished  by  further  decreases  in 
saturation  mixing  ratio  or  by  moisture  transport  into  the  cloud  as  water 
precipitates  from  the  cloud. 

The  evolution  of  the  cloud  and  precipitation  variables  for  a  range  of 
turbulence  dissipation  rates  for  the  conditions  of  marine  cloud  in  Table  A6 
are  shown  in  Figures  A5  through  A7.  The  most  basic  and  general  trend  in  these 
illustrations  is  in  the  decrease  of  the  time  to  reach  a  maximum  precipitation 
water  level  as  well  as  the  time  for  onset  of  a  significant  precipitation  flux 
as  the  turbulence  levels  rise.  This  result  is  a  manifestation  of  the 
increased  rates  N(.(.  and  N^-p  (which  control  the  auto-conversion  of  cloud)  as 
the  turbulence  dissipation  level  rises.  At  the  very  highest  turbulence  levels 
e  ~  1-10  m^/s^  the  cloud  droplets  are  progressively  de-correl ated  with  the 
smallest  scale  eddies  in  the  cascade  and  the  turbulent  shearing  mechanism 
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Precipitation  Flux  for  homogeneous  cloud  illustration  for 
various  turbulonce  dissipation  rates. 
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Figure  A. 6  Precipitation  mixing  ratio  for  homogeneous  cloud  illustration 
for  various  turbulence  dissipation  rates. 
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Precipitation  drop  radii  for  homogeneous  cloud  illustration  for 
various  turbulence  dissipation  rates. 


becomes  limited;  the  turbulent  acceleration  a  continues  to  grow,  however,  as 


TABLE  A6.  CONDITIONS  FOR  HOMOGENEOUS  CLOUD  ILLUSTRATION 


q£(0) 

1.0  g/kg 

^po 

30  )j(n 

^c 

1  km 

n^  (maritime) 

200/cm3 

It  will  be  observed  for  these  evolutions  at  fixed  initial  cloud  water 
content  that  higher  turbulence  levels  promote  higher  levels  of  maximum 
precipitation  flux  at  relatively  smal ler  average  precipitation  drop  sizes, 
this  is  because  at  fixed  water  content,  the  precipitation  drop  size  is 
controlled  by  the  ratio  of  auto-conversion  rate  N^-p  to  cloud  collection  rate 
Np^,.  At  higher  turbulence  levels,  the  turbulent  shear  contribution  to  N^-p 
enhances  this  ratio  and  hence  increases  the  number  density  of  precipitation 
drops  which  are  formed  before  the  cloud  is  completely  collected.  The  result 
of  the  enhanced  precipitation  number  density  is  a  relatively  smaller 
precipitation  drop  size  since  the  total  water  content  is  fixed. 

In  more  complex  natural  cloud  evolutions  (in  contrast  to  this  simplified 
illustration  in  which  a  fixed  amount  of  cloud  water  appears  instantaneously 
independent  of  turbulence  level)  higher  turbulence  levels  are  strongly 
correlated  with  high  updraft  levels  (or  stronger  cloud  top  radiative  cooling 
in  stratus)  and  correspondingly  higher  liquid  water  content.  Since 
precipitation  drop  size  increases  directly  with  liquid  water  content,  larger 
precipitation  drop  size  predicted  by  the  present  model  will  be  correlated  with 
higher  turbulence  levels  in  such  natural  cloud  evolutions. 

It  may  be  noted  that  at  the  highest  levels  of  turbulence  (e  -  1-10  m^/s^) 
a  significant  precipitation  flux  is  established  over  several  tens  of  minutes. 
Such  turbulence  levels  may  be  characteristic  of  the  dynamics  within  strong 
cumulus  cells  and  it  is  of  interest  to  examine  the  evolutions  predicted  here 
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with  those  of  a  popular  precipitation  model  (Kessler,  1969)  which  is 
parameterized  for  strong  cumulus  convection  and  contains  no  mechanistic 
description  of  the  role  of  turbulence.  The  results  of  the  same  case  presented 
in  Figures  A5  through  A7  predicted  by  the  Kessler  model  are  shown  in 
Figure  A8.  the  time  scales  and  general  evolution  seem  comparable  to  the 
present  model  for  turbulence  levels  £  i  Im^/sec^,  save  for  the  smaller 
precipitation  sizes  predicted.  It  should  be  noted,  however,  that  the  present 
model  has  no  empirical  parameters  restricting  it  to  such  cumulus  clouds. 
Thus,  the  high  turbulence  levels  of  the  present  model  generate  precipitation 
of  high  number  density  and  moderate  size  (Rp  ^  450  ym).  Thus,  while  the 
Kessler  model  by  virtue  of  its  parameterization  is  incapable  of  describing  the 
stratus  case  in  absence  of  updraft,  we  believe  the  present  model  when 
integrated  with  fluid  dynamic  mean  motion  including  an  updraft,  would  predict 
rain  drop  sizes  consistent  with  the  Kessler  model  and  the  Marshal -Palmer 
parameterization  which  is  the  key  part  of  its  model  structure. 

Concluding  Remarks 

Atmospheric  turbulence  can  play  a  significant  role  in  the  evolution  and 
development  of  cloud  and  precipitation.  Both  collisional  and  non-coll isional 
turbulent  processes  can  be  operative;  in  the  present  work  we  have  developed 
the  turbulent  collisional  effects  in  detail,  incorporated  them  in  a  simplified 
microphysical  closure  model,  and  carried  out  the  appropriate  turbulent 
ensemble  averages  with  attention  to  the  droplet  correlations  with  the 

turbulent  eddy  cascade.  Further  study  is  required  to  define  the  conditions 

under  which  non-col  1 isional  turbulence  mechanisms  are  important. 

The  operational  model  for  cloud  and  precipitation  which  has  been 

developed  exhibits  results  in  terms  of  time  scales  and  magnitudes  of 
precipitation  sizes  and  fluxes  which  are  consistent  with  those  naturally 

occurring  in  the  atmosphere.  The  turbulence  levels  which  effect  these  results 
are  typical  of  naturally  occurring  turbulence  levels  in  the  atmosphere. 

The  present  model  provides  an  extension  of  Kessler's  parameterization  by 
providing  a  direct  dependence  on  turbulence  level  of  his  parameter  for 
autoconversion  from  cloud  water  to  rain  water.  This  generalization  is 
accomplished  at  the  expense  of  carrying  two  more  variables;  the  turbulence 
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Figure  A. 8  Evolution  of  the  precipitation  water  mixing  ratio  for  the 
homogeneous  cloud  evolution  conditions  of  Figures  (A.5)-(A.7) 
as  predicted  by  the  Kessler  model.  There  is  no  dependence  upon 
turbulence  in  the  Kessler  Model. 


dissipation  rate  and  the  cloud  drop  number  density.  Even  if  the  turbulence 
dissipation  rate  is  not  directly  carried  in  a  numerical  cloud  model,  it  can  be 
approximately  related  to  whatever  dynamical  variables  are  provided  to  control 
entrainment.  Variations  in  the  cloud  drop  number  density  from  its 
environmental  input  value  provide  the  measure  of  dispersion  in  the  cloud 
droplet  spectrum  which  is  essential  in  describing  the  evolution  from  cloud  to 
precipitation. 
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APPENDIX  B 


INCORPORATION  OF  AN  ANISOTROPIC  LENGTH  SCALE 
INTO  SECOND-ORDER  CLOSURE  MODELING 
OF  THE  REYNOLDS  STRESS  EQUATION 

★ 

by  R.  I.  Sykes,  C.  Cerasoli  * 

W.  S.  Lewellen  and  C.  Swanson 

INTRODUCTION 

A  critical  feature  of  any  second-order  closure  model  is  how  the 
macroscopic  nature  of  a  given  turbulent  flow  field  is  incorporated  into  the 
model.  A  second-order  closure  model  attempts  to  provide  a  unique  relationship 
between  the  means,  variances  and  covariances  of  the  primary  variables 
independent  of  specific  boundary  conditions.  Since  turbulence  is  a  property 
of  the  macroscopic  flow  field  rather  than  a  local,  single  point  property  of 
the  fluid,  it  is  natural  to  expect  that  some  information  from  the  two  point 
averages  will  be  required  to  uniquely  define  the  relationships  between  the 
first  and  second-order  moments  at  a  single  point.  In  current  models  this 
macroscale  information  is  supplied  either  by  a  length  scale  equation,  or 
equivalently  an  equation  governing  the  dissipation  of  turbulent  kinetic 
energy.  All  other  macroscales  entering  the  problem  are  then  assumed  to  be 
directly  proportional  to  this  single  macroscale.  This  assumption  is  not 
universally  valid,  particularly  in  the  presence  of  a  wall  and/or 
stratification. 

Lewellen  and  Sandri  (1980)  made  an  attempt  to  introduce  a  simple  two 
scale  approach  in  order  to  permit  the  horizontal  velocity  variance  to  follow 
completely  different  scaling  than  that  obeyed  by  the  vertical  velocity 
variance  in  the  atmospheric  surface  layer  under  unstable,  convective 
conditions.  The  resulting  two  scale  model  permitted  the  horizontal  velocity 
variance  to  scale  with  the  mixed  layer  height,  consistent  with  available  data, 
rather  than  follow  the  Moni n-Obukhov  scaling  which  governs  the  vertical 
velocity  variance.  The  present  report  describes  how  this  same  result  can  be 


*The  contributions  of  C.  Cerasoli  and  C.  Swanson  to  this  work  were  supported 
by  other  Navy  contracts. 
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obtained  in  a  more  natural  way  which  provides  a  firmer  foundation  for 
extending  the  model  to  more  general  flows.  By  considering  the  turbulent 
energy  to  be  composed  of  two  populations  of  eddies,  one  large  and  one  small, 
the  two  length  scales  appear  naturally.  The  only  empirical  information  needed 
to  complete  the  model  is  an  algebraic  relationship  governing  the  partition  of 
the  Reynolds  stress  between  the  large  and  the  small  scales. 


A  MODEL  FOR  THE  WALL-EFFECT 


Our  conceptual  model  of  the  turbulent  flow  near  a  wall  without  shear  is 
motivated  by  the  turbulence  spectra  measured  near  walls.  The  spectra  of 
Thomas  and  Hancock  (1977),  Willis  and  Deardorff  (1974),  and  Kaimal  (1978)  all 
show  the  same  general  features,  namely  that  the  wavelength  of  the  peak  in  the 
normal  velocity  component  spectrum  decreases  linearly  to  zero  at  the  wall, 
while  that  of  the  tangential  components  remains  roughly  constant.  Kaimal 's 
spectra  are  shown  in  Figure  Bl.  This  is  interpreted  straightforwardly  as 
large  eddies  being  forced  to  flow  tangential  to  the  wall,  so  that  their  normal 
velocity  component  vanishes,  and  only  small  local  eddies  contain  any 
significant  normal  energy. 

Previous  second-order  closure  models  only  contain  one  length  (or  time) 
scale,  and  this  always  goes  linearly  to  zero  at  the  wall.  Thus  we  expect  such 
models  to  be  able  to  predict  normal  component  correlations  reasonably  well, 
but  not  the  horizontal  ones.  Our  remedy  for  this  problem  is  to  consider  the 
flow  near  the  wall  to  be  comprised  of  two  populations  of  eddies.  The  first  is 
the  small  scale  population  which  is  the  one  currently  described  by 
second-order  closure  models.  This  population  has  a  length  scale  proportional 
to  the  distance  from  the  wall,  and  is  fully  three-dimensional,  i.e.,  all 
components  of  the  Reynolds  stress  tensor  have  the  same  order  of  magnitude. 
The  second  population  are  the  large  eddies  with  a  length  scale  determined  by 
the  flow  scales  away  from  the  wall.  These  eddies  are  two-dimensional  at  the 
wall  since  their  normal  energy  component  vanishes.  We  therefore  write  our 
modeled  Reynolds  stress  equations  as  follows: 
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Figure  Bl.  Normalized  velocity  spectra  as  a  function  of  height  for  the 
atmospheric  boundary  layer  under  themally  unstable  conditions 
(Kaimal ,  1978).  The  spectral  peaks  for  the  horizontal 
velocities,  u  and  v,  remain  constant,  while  the  spectra  peak 
for  w  moves  to  small  lengths  as  the  ground  is  approached. 
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where  u^Uj  is  the  total  Reynolds  stress,  and  u^uj^ ,  u^Uj^  represent  the  large 
and  small  contributions,  respectively;  also  qf=u-u,-*-,  and  q2=u  -u.^.  A,  is 
the  normal  or  vertical  (anticipating  the  convective  application)  length  scale, 
and  is  the  tangential  or  horizontal  length  scale.  u^0  is  the  turbulent 
heat  flux,  g^-  is  the  gravitational  vector,  and  n^-  =  (0,0,l)  is  the  unit  vector 
normal  to  the  wall.  Thus  our  model  consists  of  a  three-dimensional  "return  to 
isotropy"  term  for  the  small  eddies,  and  a  corresponding  two-dimensional  term 
for  the  large  eddies.  The  dissipation  term  depends  only  on  the  small  eddies, 
since  the  large  eddies  must  cascade  their  energy  through  the  small  scale 
eddies  before  it  is  ultimately  dissipated.  The  diffusion  term  contains  small 
eddies  diffusing  isotropically,  whilst  large  eddies  only  diffuse  in  the 
tangential  plane.  Note  that  if  u-jUj^=0,  i.e.,  all  the  energy  is  in  small 
scales,  then  we  recover  the  single  scale  equations  of  Lewellen  (1977)  with 
b=l/8,  V(,=0.3. 

At  this  stage  we  have  not  specified  the  partition  of  u^uj  between  large 
and  small  scales;  we  shall  make  this  partition  algebraically  by  reference  to 
the  measured  spectra.  The  horizontal  energy  spectra  from  the  experiment  of 
Thomas  and  Hancock  (1977),  show  that  after  the  energy  peak  at  a  scale 
determined  by  the  free  stream  turbulence,  the  spectrum  falls  off  smoothly  into 
an  inertial  range  with  no  irregularity  at  the  smaller  scales  where  the 
vertical  energy  peaks.  This  spectral  shape  is  very  similar  in  Kaimal's 
atmospheric  measurements--provided  the  spectra  are  taken  outside  the  surface 
shear  layer,  i.e.,  z/L  ~  5,  where  z  is  the  height  above  the  ground  and  L  is 
the  Monin-Obukhov  length.  The  existence  of  an  inertial  range  prompts  the  use 
of  a  simple  scaling  law  to  estimate  the  small  scale  horizontal  energy,  viz. 
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Thus,  the  horizontal  energy  at  wavelengths  of  order  are  simply  scaled 
down  from  the  large  wavelength  total  by  the  ratio  of  length  scales  to  the 
two-thirds  power.  To  complete  the  model,  we  have 


UiUj  ^i^j> 


i=3  or  j=3 


and 


Thus  al  1  the  energy  in  correlations  involving  the  vertical  component  is 
contained  in  small  scales,  and  the  large  eddy  energy  is  the  remainder;  hence 
there  is  no  large  eddy  energy  in  vertical  components. 

To  complete  the  model  specification,  we  write  the  heat  flux  equations  as 


D 

Dt 


“77  -  9T 

u^-O  =  -  u.u^-  - — 
1  1  J  9x,- 


(B.2) 


and  temperature  variance: 


^  02  =  -  2ui0  -  2bs  ^  02  (B.3) 
Dt  J  8xj  Ay 

These  are  the  standard  equations  from  Lewellen  (1977),  and  we  use  the 
same  coefficients,  i.e.,  A=0.75,  s=1.8,  which  is  tantamount  to  assuming  that 
the  temperature  fluctuations  are  all  in  the  small  scale  part  of  the  spectrum. 
This  is  justified  because  the  production  terms  in  these  equations  only  involve 
the  vertical  velocity  component  which  is  only  present  at  small  scale. 

Finally,  length  scales  are  specified  by  the  same  scale  equation  as 
Lewellen  (1977)  for  the  vertical  scale,  but  using  small  scale  energies 
appropriately,  i.e. 
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and  =  max  where  the  maximum  is  taken  over  the  domain. 


(B.4) 


Thus,  in  essence  we  accept  that  the  original  model  described  by 
Lewellen  (1977)  does  a  creditable  job  of  predicting  the  small  scale  part  of 
the  spectrum  near  the  wall,  and  we  therefore  maintain  the  small  scale 
dynamics.  However,  we  include  a  two-dimensional  large  scale  component  which 
is  dominant  near  the  wall,  and  use  an  inertial  range  scaling  law  to  obtain  the 
partition  between  the  two  energies.  Note  that  although  this  conceptual  model 
was  arrived  at  through  consideration  of  the  turbulence  close  to  the  wall,  the 
model  reverts  naturally  to  the  single  scale  model  as  one  moves  out  to  the 
middle  of  the  flow  because  A^  -  A^  in  the  central  region,  and  therefore 


MODEL  RESULTS  FOR  SHEAR-FREE  FLOWS 


Convection  between  Flat  Plates  at  Fixed  Temperature 

The  first  experimental  flow  we  compare  with  is  thermal  convection  between 
horizontal  flat  plates  held  at  fixed  temperatures.  The  data  for  comparison 
are  from  Deardorff  and  Willis  (1967),  and  we  use  their  highest  Raleigh  number 

7 

data  of  10  ,  since  we  include  no  low  Reynolds  number  terms  in  our  model 
equations.  Equations  B.l  through  B.4  were  solved  numerically  with  an  imposed 
heat  flux  at  the  lower  boundary,  z=0,  and  plane  of  symmetry  conditions  at 
z=D/2  where  D  is  the  distance  between  the  horizontal  plates.  At  the  lower 
boundary  we  set 


Bo  ^  ^  r\9  n 

-  U^  =  -  V^  =  -  0^  =  0 

9Z  8Z  3Z 


and  w2  =  0,  w0  =  Hg. 
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Figure  B2  shows  the  r.m.s.  turbulence  components  from  the  model  scaled  by 
w*  where 


wl  =  ^  HqD. 

T  ° 

The  experimental  values  are  plotted  also  as  continuous  lines  with  error  bars 
indicating  the  general  level  of  scatter  in  the  data.  Note  that  the 
experimental  values  had  to  be  rescaled  for  this  plot,  and  it  was  necessary  to 
use  a  definite  value  for  the  heat  flux  to  achieve  this.  The  value  used  was 
the  Globe  and  Dropkin  (1959)  empirical  function  value  since  Deardorff  and 
Willis  regard  their  direct  measurements  of  turbulent  heat  flux  as  less 
reliable;  the  difference  between  the  values  is  somewhat  less  than  10%.  There 
is  good  agreement  between  the  predicted  and  measured  values  for  both  vertical 
and  horizontal  components.  There  is  a  slight  increase  in  horizontal  energy 
near  the  wall  in  the  experiment,  and  the  magnitude  of  this  increase  is 
underpredicted  but  the  model  does  maintain  the  horizontal  energy  right  up  to 
the  wall  and  does  actually  show  a  10%  increase  from  the  value  in  the  middle  of 
the  flow.  Below  z=0.075D,  the  measured  energy  falls  off  toward  zero,  which  is 
an  indication  of  the  viscous  boundary  layer  effect.  The  vertical  energy  level 
is  also  well  predicted,  with  some  over  prediction  for  small  z,  which  again  is 
probably  a  viscous  effect. 

Figure  B3  shows  r.m.s.  temperature  fluctuations  scale  by  0*=Hq/w*,  and 
again  the  predicted  values  are  very  accurate  except  in  the  region  close  to  the 
surface  where  molecular  diffusion  becomes  significant. 

Convection  from  a  Heated  Plate  with  an  Overlying  Stable  Layer 

In  this  flow,  a  constant  heat  flux  Hq  is  applied  at  the  base  of  an 
infinite  layer  of  stably  stratified  fluid,  and  a  convective  mixed  layer  grows 
in  depth  with  time.  The  model  was  run  until  a  self  similar  solution  was 
obtained,  and  the  results  were  plotted  in  dimensionless  form  with  z^  as  the 
length  scale,  and  w4  =  gHQZ^-/T  as  velocity  scale.  Here  z^-  is  the  depth  of  the 
mixed  layer,  measured  to  the  point  with  maximum  temperature  variance  at  the 
top  of  the  layer.  The  experimental  values  are  from  Willis  and 
Deardorff  (1974). 
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Figure  B2.  R.m.s.  turbulent  velocity  coiiiponents  for  corivection 
flat  plates. 

Model  predictions;  —  horizontal  component 

—  vertical  component 

Laboratory  data  from  Deardorff  and  Willis  (1967); 

horizontal  component 
vertical  component 

Error  bars  indicate  the  scatter  in  the  data. 


between 
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Figure  B3.  R.m.s.  temperature  fluctuations  for  convection  between  flat 
plates.  Solid  line  represents  the  laboratory  data,  and  the 
model  predictions  are  plotted  as  dots.  Error  bars  indicate  the 
scatter  in  the  data. 
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Figure  B4  shows  velocity  variances,  and  excellent  agreement  between 
measured  and  predicted  profiles.  It  is  of  interest  to  note  that  in  this  case 
the  horizontal  variance  is  slightly  over  predicted  at  the  surface.  Figure  B5 
shows  the  temperature  variance  which  is  very  closely  predicted  in  the  lower 
80%  of  the  mixed  layer,  but  the  large  maximum  at  the  inversion  is  not 
predicted.  This  is  a  failure  of  the  second-order  closure  model  to  predict  the 
entrainment  fluxes,  which  Zeman  and  Lumley  (1976)  have  noted,  but  is  not  the 
subject  of  this  paper. 

Shear-free  Turbulent  Boundary  Layer 

In  this  experiment  by  Thomas  and  Hancock  (1977),  a  turbulent  free  stream 
is  made  to  pass  a  plane  boundary  moving  at  the  same  mean  speed.  This 
experiment  is  a  rather  inconclusive  test  of  the  model  for  two  reasons. 
Firstly,  the  primary  effect  of  the  moving  wall  is  the  generation  of  an 
impulsive  pressure  field  which  brings  the  normal  velocity  to  rest  at  the 
boundary,  and  distorts  the  turbulence  spectra  as  described  by  Hunt  and 
Graham  (1978).  Our  second-order  closure  model  cannot  predict  this 
instantaneous  change,  and  therefore  we  need  to  initialize  the  calculation  with 
a  field  that  satisfies  the  boundary  conditions.  The  experimental  section  was 
not  very  long,  so  the  results  are  dependent  on  the  initial  conditions  and  the 
latter  were  not  measured.  Secondly,  the  experiment  is  clearly  non-ideal  in 
the  sense  that  the  two  tangential  components,  streamwise  and  transverse, 
behave  quite  differently  with  downstream  distance;  the  two  components  must  be 
identical  in  the  ideal  experiment,  since  there  is  no  mean  flow  in  the  frame  of 
the  moving  wall.  Furthermore,  we  found  the  experimenter's  suggestion  that 
streamwise  pressure  gradients  were  responsible  for  the  differences  to  be 
unsubstantiated,  since  these  were  of  much  too  small  a  magnitude  to  effect  any 
significant  changes  when  included  in  the  calculation. 

We  initialized  the  model  with  the  theoretical  profiles  of  Hunt  and 
Graham,  and  then  integrated  forward  approximately  six  turbulence  time  scales, 
A/q,  where  A  and  q  refer  to  free  stream  values.  This  integration  distance  was 
estimated  from  the  measurements  of  turbulent  energy  and  length  scales  at  the 
downstream  portion  of  x/M=25.  In  this  comparison  we  have  used  the  relation 
Ly=1.5A  to  relate  the  integral  length  scale  of  the  measurements  to  the  model 
length  scale.  This  relation  was  obtained  by  Sandri  (1977)  for  isotropic 
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Figure  B4.  Turbulent  velocity  variances  for  inversion-capped  convective 
layer.  Symbols  are  as  in  Figure  Bl,  and  the  laboratory  data  is 
from  Willis  and  Deardorff  (1976). 
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Figure  B5.  Temperature  variances  for  inversion-capped  convective  layer. 
Symbols  as  in  Figure  B2. 
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turbulence  and  is  valid  as  we  use  it  only  for  the  free  stream  values. 

Figure  B6  shows  measured  and  predicted  profiles  of  the  normal  energy, 
and  the  transverse  energy,  v^.  Velocities  are  scaled  by  the  free  stream 
turbulent  velocities,  and  the  normal  coordinate  is  scaled  by  the  free  stream 
integral  length  scale.  We  note  that  in  the  experiment,  shows  a  strong  peak 
at  the  wall  which  increases  with  downstream  distance.  Hunt  and  Graham  suggest 
this  could  be  some  disturbance  close  to  the  wall  diffusing  outward,  and  our 
predictions  undoubtedly  show  much  better  agreement  with  the  transverse 
component,  since  u^=?v2  in  the  model.  The  main  features  of  the  prediction  are 
that  the  initial  peak  of  1.5  at  the  wall  in  v^  is  rapidly  reduced,  and  has 
fallen  to  1.03  at  x/M=25.  The  asymptotic  model  value  is  0.95  far  downstream. 
The  prediction  of  v^  is  in  close  agreement  with  the  measurements,  which  show  a 
virtually  constant  value.  However,  the  model  prediction  of  the  development  of 
W^  is  that  the  initial  profile  diffuses  oufward,  while  the  measurements  show 
an  pned^nging  profile*  Thus  at  x/M=25,  we  have  the  correct  shape  for  the 
profile,  but  the  scaling  length  appears  to  be  too  big  by  almost  a  factor  of  2. 
In  view  of  the  nop-ideal  nature  of  the  experiment  in  terms  of  the  behavior  of 
u^,  it  seems  unwise  to  speculate  on  the  discrepancies  in  the  comparisons.  All 
we  can  really  conclude  is  th^t  the  model  does  predict  the  observed  qualitative 
results. 


THE  EFFECT  OF  WALL-LAYER  SHEAR 

In  the  atmospheric  boundary  layer,  we  generally  have  both  thermal 
convection  and  mean  shear  present  at  the  same  time.  The  relative  importance 
of  these  two  effects  in  the  main  part  of  the  boundary  layer  is  measured  by  the 
rati  0  L/z-j ,  where 

u^T 

L  =  — —  ,  the  Monin-Obukhov  length, 

<gHo 

and  z^  is  the  inversion  height;  u*  is  the  surface  friction  velocity,  Hq  the 
surface  heat  flux,  and  k  is  von  Karman's  constant,  taken  to  be  0.39  here. 
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Figure  B6.  Turbulent  velocity  variances  for  the  shear-free  wall  layer 
compared  with  the  data  of  Thomas  and  Hancock  (1977).  Symbols 
as  in  Figure  Bl. 
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When  z  <  0(L),  shear  effects  become  important  and  our  simple  model  of  the 
partition  between  large  and  small  eddies  is  then  no  longer  appropriate.  In 
this  region,  the  shear  production  term  in  the  turbulence  energy  equation  is 
directly  producing  significant  small  scale  horizontal  energy,  and  our  inertial 
range  scaling  law  will  fail.  This  can  be  seen  in  the  spectra  of  Kaimal 
(1978),  which  show  a  bulge  at  higher  wave  numbers  in  the  horizontal  spectra 
where  z  <  0(L).  We  therefore  require  a  more  sophisticated  energy  partition. 
We  have  tried  to  keep  the  partition  as  simple  as  possible,  since  our  principal 
goal  in  this  paper  is  the  development  of  a  model  for  large  eddies  near  a  wall. 
We  shall  describe  an  extension  of  the  algebraic  partition  which  allows  the 
computation  of  shear  flows,  and  although  we  base  the  model  on  a  crude  physical 
model  of  the  energy  spectrum,  it  should  be  regarded  mainly  as  a  simple  device 
to  extend  the  range  of  applicability  of  the  model  to  cover  practical 
calculations  of  the  atmospheric  boundary  layer. 

The  alternative  to  an  algebraic  partition  would  be  to  carry  two  sets  of 
equations  for  large  and  small  scales  separately,  in  the  manner  introduced  by 
Hanjalic  et  al .  (1979).  Ultimately  this  may  well  be  the  best  way  to  proceed, 
but  the  interaction  terms  between  the  two  scales  require  a  good  deal  of 
modeling  effort,  because  Hanjalic  et  al.'s  assumption  of  isotropy  in  the  small 
scales  is  inappropriate  here.  We  therefore  prefer  to  keep  the  model  as  simply 
as  possible  at  this  stage,  and  at  least  determine  what  conditions  a  more 
general  model  needs  to  satisfy. 

Our  assumption  will  be  that  the  energy  in  the  small  scales  arises  both 
from  the  cascade  from  larger  scales  and  also  from  direct  production  in  the 
small  scale  part  of  the  spectrum.  The  small  scale  production  in  the  case  of 
the  atmospheric  boundary  layer  is 


9U: 


^s  "  ^i^j  8x/ 


i.e.  the  shear  production  term,  since  this  is  generating  small  scale  energy  in 
the  horizontal  components  near  the  wall.  Since  this  becomes  increasingly 
dominant  near  the  wall  where  the  time  scales  are  very  short,  we  have 
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i.e.  the  production  balances  dissipation.  In  fact,  these  terms  are  so 
singular  in  the  logarithmic  layer,  that  it  is  necessary  to  allow  this  exact 
balance  very  close  to  the  wall.  Thus 
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(B.5) 


Before  combining  this  result  with  the  energy  cascade  result,  a  further 
complication  of  the  shear  production  is  that  it  does  not  necessarily  generate 
u|  and  v|  equally.  A  simple  equilibrium  argument  shows  that 


1  +  Psu/Ps 


u|  +  v| 


where  is  the  shear  production  term  in  the  equation.  Thus,  defining 


1  "  Psu/Ps 


and 
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so 
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our  final  model  for  the  small  scale  horizontal  components  is 


u|  =  a  fgQ  (u^+v^)  + 


73 


so 


>  0 


(B.6) 


v|  =  (1-a)  fgo  (u2+v2)  + 


(B.7) 


Otherwi se 


(B.8) 


Thus  we  have  a  simple  algebraic  partition  which  satisfies  the 
requirements  of  allowing  an  exact  small  scale  balance  very  close  to  the  wall, 
and  going  back  to  our  original  inertial  large  scaling  when  there  is  negligible 
small  scale  production. 

RESULTS  FOR  THE  ATMOSPHERIC  SURFACE  LAYER 


The  main  point  of  our  large  eddy  model  is  that  the  energy  in  the 
horizontal  components  near  the  surface  is  "passive",  in  the  sense  that  it  does 
not  affect  the  vertical  correlations.  This  is  indeed  the  case,  and  all  the 
similarity  variables  in  the  surface  layer  are  within  a  few  percent  of  the 
results  of  Lewellen  and  Teske  (1974)  and  are  therefore  not  presented  again 
here.  However,  the  horizontal  components  now  show  a  larger  energy  depending 
on  the  value  of  the  horizontal  scale,  A^.  In  the  following  comparisons 
A|^=0.235z^;  this  result  is  obtained  from  the  scale  Equation  B.4. 


Figure  B7  shows  a  typical  surface  layer  profile  for  the  three  energy 
components  with  z.j/L=50.  The  velocities  are  scaled  by  the  surface  friction 
velocity,  u*.  The  features  to  note  are  that  w^  is  1.4  at  the  wall,  in  accord 

2  /o 

with  the  neutral  log  layer  value,  and  then  increases  asymptotically  like  z 
for  large  z.  On  the  other  hand,  the  horizontal  components  are  larger,  and 
remain  relatively  constant  right  down  to  the  wall.  There  are  variations  of  a 


161 


few  percent  close  to  the  wall,  and  this  is  due  to  the  precise  nature  of  our 
small  scale  energy  partition  function.  The  atmospheric  measurements  show  the 
horizontal  energy  remaining  constant  within  the  scatter  of  the  data  which  is 
about  10%.  We  might  note  that  several  simpler  versions  of  the  partition 
function  (i.e.,  not  satisfying  all  the  requirements  of  Section  4)  all  gave 
unreasonable  profiles  in  the  sense  of  energy  doubling  at  the  wall  or 
decreasing  by  50%. 

A  comparison  of  the  r.m.s.  horizontal  energy. 


against  the  observed  values  is  shown  in  Figure  B8.  We  have  plotted  the  value 
of  from  the  model  at  z/L=3,  i.e.,  outside  the  region  of  the  surface 

variation.  The  model  predictions  lie  well  within  the  scatter  of  the 
observations  at  moderate  z^-/L;  there  are  only  two  data  points  at  high  values, 
and  the  model  over  predicts  these  by  about  10%. 

Finally,  we  note  that  with  the  formulation  of  small  scale  energy  given  by 
Eqs.  B.5-B.7,  al  1  the  energy  is  contained  in  the  small  scales  for  neutral 
flows,  e.g.,  boundary  layer  or  channel  flows,  and  also  in  the  stable 
atmospheric  surface  layer.  Hence,  for  these  flows  the  predictions  are 
unchanged  from  the  earlier  single  scale  model. 


162 


3 


N 

5 

eg 

> 


M 

3 


163 


Figure  B7.  Predicted  turbulent  velocity  variances  in  the  atmospheric 
surface  layer  for  z.j/L=50. 
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Figure  B8.  Comparison  between  r.m.s.  horizontal  turbulence  energy 
predictions  and  observations  as  a  function  of  z-j/L.  T!ie  model 
prediction  is  shown  as  a  solid  line.  The  data  aie  from 
Panofsky  et  al.  (1977). 


CONCLUSION 


A  second-order  closure  model  which  accounts  for  the  effects  of  a  rigid 
wall  on  large  scale  turbulent  eddies  has  been  presented.  The  basis  of  the 
model  is  a  partitioning  of  the  turbulent  kinetic  energy  between  small  scale 
three-dimensional  eddies  and  large  scale  eddies  which  only  have  tangential 
energy  components  near  the  wall.  The  relative  partitioning  of  the  tangential 
energy  has  been  made  using  the  inertial  range  spectral  distribution.  The 
model  performs  well  for  a  number  of  flows  in  which  the  energy  is  principally 
in  large  scale  eddies. 

A  somewhat  over-simplified  description  of  the  model  would  be  that  we 
retain  the  original  single  scale  equations  for  the  small  scale  component,  and 
then  scale  up  tangential  components  using  the  ratio  of  the  two  length  scales. 
This  view  has  some  utility  in  that  it  is  correct  except  for  the  turbulent 
diffusion  terms  in  the  Reynolds  stresses;  it  also  has  some  basis  in 
observation,  since  the  spectra  of  tangential  velocity  components  do  show  an 
apparently  undisturbed  inertial  range  fall  off  from  the  energy  containing 
length  scales.  However,  the  diffusion  terms  are  not  always  negligible  in  the 
flows  we  have  considered,  so  that  there  is  some  effective  coupling  between  the 
two  scales;  in  fact  the  diffusion  is  the  dominant  mechanism  in  the  surface 
layer  profiles  of  horizontal  energy  in  Section  4.  Furthermore,  the  dynamical 
basis  of  this  model  is  very  useful  in  developing  corresponding  two  scale 
equations  for  passive  scalars  in  the  atmosphere,  where  the  small  scale  is  no 
longer  the  controlling  scale  in  the  equations;  the  work  on  scalar  diffusion 
will  be  presented  in  a  subsequent  paper. 

The  main  conclusion  from  this  study  is  that  several  aspects  of  wall 
turbulence  can  be  accurately  described  using  the  very  simple  conceptual  notion 
of  passive  horizontal  energy,  so  that  the  active  turbulence  is  calculated  with 
the  usual  closure  model,  but  the  active  part  of  the  turbulence  is  derived  from 
the  total  using  a  simple  hypothesis  for  the  spectral  shape. 
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